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化 学 工程 师 的 任务 
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T UL 








化 学 工程 专业 数学 模型 类 型 
非 线 性 方程 (组 ) 


Te 


第 微分 方程 (组 ) 解析 解 ， 必须 采用 数值 
解法 
模型 的 数值 解法 是 应 用 
数学 的 一 个 分 支 ， 通 第 
偏 微分 方程 (组 ) 称 为 计算 数学 (数值 分 
术 ， 数 值 方法 ) 











化 学 工程 第 用 软件 


数学 软件 : 
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WELLecte 
WE 


ieSiezl 


化 工 模 拟 软件 : 
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盖 AspenPlus 
=] IT@7A DB 
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本 课程 的 学 习 目 的 


对 于 数值 分 析 
的 内 容 个 过 多 
涉及 ， 只 注意 
数 信 计算 结 朱 
的 准确 性 


尝 会 Matlab 的 使 用 ， 可 以 利 
用 Matlab 求 解 较为 复杂 的 化 
工 数学 模型 








化 工 专业 知识 
作为 背景 ， 不 
涉及 模型 的 推 
导 ， 注 重 模型 
求解 过 程 的 方 





本 课程 基本 内 容 


第 一 讲 Matlab 简 介 与 基本 数学 运算 

第 二 讲 非 线性 方程 组 求解 与 迭代 法 
第 三 讲 天 阵 操作 与 线性 方程 组 求解 

第 四 讲 插值 、 拟 合 与 数值 微分 、 积 分 
第 五 讲 常 微 分 方程 数值 解 

第 六 讲 偏 微分 方程 数值 解 

第 七 讲 统计 初步 与 最 优化 方法 





学 习 本 课程 的 注意 事项 吧 





Ge- 


*。 学 好 本 课程 的 唯一 途径 是 王 上 上 人 记 实 跨 


。 数值 计算 效率 和 效果 的 保证 有 很 多 技巧 ,可 以 参考 数 
从 方法 【六 从 分 析 ) 方 的 和 


DO 刘 则 壮 ， 本 技术 与 Matlab， 科 学 出 版 社 


[| 同 洲 大 学 计算 骂 缉 研 窜 ， 现代 数值 数学 和 计 
算 ， 村 祁 大 这 引 反 往 


D 黄 华 江 ， 实 用 化 工 计 算 机 模拟 ， 化 学 工业 出 版 社 


张 志 涌 ， 精 通 Matlab6.5 版 ， 北 京 航空 航天 大 学 出 
版 社 


。 对 于 数值 计算 的 结果 ， 应 注意 分 析 结 果 的 意义 





ELE|eEE 轴 CO 


Matlab 喜 Wairix Lapotary 的 
缩写 ， 最 初 是 美国 新 墨 西 可 
大 学 Moler 教 授 编写 的 
LINPACK 和 EISPACK 接 口 
程 订 





1984 年 ，MathWorks 公 司 创 
建 ，MATLAB 正 式 推 向 市 场 


20 世 纪 90 年 代 以 来 ，MATLAB 已 
成 为 数值 计算 软件 的 佼佼 者 





9 全 [Ne [at Jack Little 











开始 的 问题 


计算 在 1/2 英 寸 不 锈 钢管 中 ， 以 2000Ib/hr 流 量 输送 水 ， 当 水 的 
温度 为 10、20、30、40、50、60、70、80C 时 ， 压 降 分 别 为 
多 少 ? 


牛 撕 流 体 在 不 锈 钢管 中 的 流动 压 降 可 由 下 式 估算 : 本 下 全 全 


其 中 ， 摩 探 压 降 ，psiy(100 英 尺 等 量 管 长 );) M， 质 量 流量 ，lb/hr; Qu ， 煌 
度 ，cP; p ， 密 度 ，Ib/ft3，D， 管 径 ，inch。 


流体 密度 可 由 下 式 档 术 :2 


p ，g/ml; 对 于 水 ，A= 0.34710; B= 0.2740; Tc=647.13K; n= 0.28571. 


TD 


HU，CcP; 对 于 水 ，A=-10.2158; B=1.7925E3; C = 1.7730E-2; D=-1.2631E-05. 








-上 本 ATIAE 
Rile 了 dit Debug 了 Desktop 业 ndow Help 


口 脓 |% 恩 家 忆 心 | 功 晤 四 |? 


















C:AProgram 了 FilesAIATLABT1Awork 用 本 
Shorteuts 加 How to k&dd 回 What s New 四 WGFforLOTFDT 
本 a Xj||  X 
| 省 男 齐 咏 久 | 渔 | 国 - 
Tane | ras | class To get started，select 了 了 TITLaB Help or Demos from the Help menu。 
>> 
由 orkspace | Currespt Directory 
NT 全 区 
当前 路 径 
和 














Matlab 的 通用 命令 


命令 说 明 | 命令 | 命令 说 明 
一 一 显示 或 改变 工作 目录 | dir “| 显示 目录 文件 
显示 文件 内 容 清除 内 存 变量 


clf 清除 图 形 窗口 pack | 收集 内 存 碎片 ， 扩 大 
内 存 空间 
cl 清除 命令 窗口 内 容 命令 窗口 信息 显示 开 




















辣 开 保持 天 四 


显示 搜索 目录 保存 内 存 变量 到 指定 
文件 
加 载 指定 文件 变量 日 志文 件 命令 
AS 


YYT 


whos | 变 有 看 | | 














通过 Help 学 习 Matlab 





在 命令 窗口 中 键入 >> help， 则 显示 以 下 内 容 : 


ULEILEleXe[s ist EL ii ellieieriewelilliitctleiS 

EleXejeiS @jeisjiilelkEUleESes wwiltUwreltS 
ELLELeNTie elelgc1lliiliie[ctiieglctelswelitillg[wt 
matlab\elmat IEUNULLEUULwSEULedLtUUEUTe IEEILe 
ELEleNa ii LECINULELUEITLILeIeIiES 

ELELeXS eeill Sellieicl[px=egliitctlmigLieteIiiS 
EleMUEIULIUI - Matrix functions - numerical linear algebra. 
ULEiLEleNeticu - Data analysis and Fourier transforms. 
ULLleNeieli taleiejicitleEcUliegeelNieulici 
ULEiLEiLeNi eeliEiylitei[eliL 汪 LeE@1B1EE:ielNsl ES 
ULEILELeNS ec - Sparse matrices. 

ULleXSeiles TITeltctileliEtiiea alte dl 二 elillire 
ULELLleXNeltclell2e 国 RLeeliiIs ileliiEEeLtclelitS 


ELcleXeLEielikie elLsiatSi[elitiEeltcTeliS 


Help + 主题 名 称 


>> help ops 
@jel:icttelEEUeESjeisielcIEeiTcUetsits 








LULUTIeeEeleisltctteltS 


plus - Plus 十 
Uplus - Unary plus 十 
minus - Minus 
uUminus - Unary minus 
mtimes - Matrix multiply 四 
times - Array multiply 
则 ee - Matrix power 和 
le 人 Lee 
ULelielz - Backslash or Ileft matrix divide \ 
Lieliei> SEEIEOITelil NELL 人 elNie[ / 
[elNiei: 区 本 21082NAelie[ 人 


elNde[ eVNAelNdeLs 岂 


数组 元 除 





-| 


Help + 函数 名 





6 必 eIeN =] 
Eee 


和 =X.^Y denotes element-by-element powers. 入 
and Y 


must have the same dlImenslons Unless one IlS a 
SCcalar. 


eic1EesUEeljeisiicitcalilteEUALILITeE 


C = POWERI(A,B) is called forthe syntax 'A .^B' 
whenAorBis an 


elej [=Ien 


Help+ 函 数 名 可 获得 详细 的 函数 使 用 方法 


Matlab 语 言 的 标点 


Ge 





珊 直 


问 量 和 和 矩阵 的 多 小 数 点 及 结构 体 
种 功能 抠 的 访问 


区 分 行 及 取消 行 续 行 符 
显示 
区 分 列 及 函数 参 注释 符 ， 百 分 号 
数 分 隔 符 
人) 指定 运算 过 程 的 Ca 


人 


加 构成 单元 数组 字符 串 标 示 符 





.分 类 方法 一 


(系统 默认 类 型 ) 


- 双 精 上 度 型 
- 单 精 度 型 
- 带 符 号 整数 
-无 符号 整数 


.分 类 方法 


ee 


三 
2 


示 
_ 向 量 


一 


-数组 


.分 类 方法 


一 A 
ES 
SS 


-实数 
-复数 





数值 的 表示 
以 下 表 达 方 式 均 合法 : 


345 -99 0.01 
1.3e-3 4.D633 
2 [1;2;3] 


[1 2; 2 11] 
3+3 6-8| 


计算 以 下 表 过 却 的 信 : 
[1 23]* [3 2 1] 








Matlab 的 计算 器 功能 


>> 2000^1.8*(10^(-10.2158+1.7925e3/283+1.773e-2*283- 
1.2631e- 5*283^2))^0.2/(20000*0.5^4.8*(0.3471*0.274A(-(1- 
283/647.13)^0.28574))/0.2323) 


加 和 可 以 得 到 绽 


ans = 
2871.8245 


命令 的 窗口 的 快捷 键 和 


快捷 键 作用 作用 

| ，Crt+N 回调 下 一 行 移 至 行 首 
一 ，Crtt+B “| 回 移 上 一 字符 移 至 行 末 
一 ，CrtttF “| 前 移 下 一 字符 删除 一 行 


Crt+ 一 左 移 一 单词 从 光标 删除 至 行 
下 


Ctrl+C 终止 正在 运行 的 
程序 








数学 函数 (elfun ) 


ET 
7 











计算 以 下 表达 陈 的 值 : 


C 〇 
吕 
吕 
Ra 
| 
一 
en 
口 一 
吕 局 工 司 9 
所 互 豆 区 页 页 
二 .二 XSDO 
六 DOS5S 
一 人 他 于 万 感 


对 
他 


format 和 人 





[ET 
[CC 
[CC 





程序 的 组 成 





数据 输入 一 运算 一 一 。 数据 输出 
键 得 输入 数学 运算 屏 姑 输出 
文件 输入 关系 运算 文件 输出 
迪 辑 运算 图 形 答 出 


流程 控制 


7 
ha 





。 芝 量 的 命名 方式 : 
一 葵 量 名 由 字母 、 数 字 和 下 划 线 组 成 ; 
一 杰 量 名 中 的 更 文学 母 六 小 号 是 有 区 列 的 ; 
一 杰 量 名 的 最 大 长 度 是 有 规定 的 
。 不同 厂 本 的 系统 规定 不 同 : 19 个 字符 、31 或 63 个 
宇 符 等 
。 可 调用 namelengthmax 函 数 得 到 系统 规定 长 度 




















变量 的 使 用 


>>clear % 删 除 工 作 区 中 所 有 定义 过 的 变量 
>>wWhos 9% 和 查看 当前 工作 区 内 变量 信息 ， 无 显示 表示 没 
有 定义 的 变量 
>> XyY=1; YX=2; % 对 变量 赋值 
>> XVy % 碍 看 变量 xy 的 当前 数值 
XV = 
1 

>> Whos 

Name Size Bytes Class 

XVy 1X1 8 double array 

VYX 1X1 8 double array 
Grand total is 2 elements using 16 bytes 
>> Clear Xy YX % 删 除 变 量 Xxy 及 yxX 
>> Whos 
>> XVy % 这 时 变量 xy 已 经 不 存在 了 


?3?3? Undefined function or variable 'Xy . 


MATLAB 系 统 的 特殊 变量 和 常数 


尹 芭 > 量 
| ans | 如 果 用 户 未 定义 变量 名 ， 系 统 用 于 计算 结果 存储 的 默认 变量 名 
pi | 圆周 率 r (= 3.1415926.…) 
无 穷 大 值 
浮 点 运算 的 相对 精度 2^(-52) 


eps 
_ realmax | 最 大 的 正 浮 点 数 ，2^((020-1 





realmin | 最 小 的 正 浮 点 数 ， 2^C1022) 
nargin | 函数 输入 参数 人 数 
nargout | 函数 输出 参数 人 数 
Iasterr | 存放 最 新 的 错误 信息 
Iastwarn | 存放 最 新 的 管 告 信息 





MATLAB 数 据 类 型 





。 数值 〈 标 量 ， 回 量 ， 数 组 ) 
。 了 字符 串 

。 单 元 数组 〈cell array ) 

。 结 构 体 〈structure ) 

。 国 数 句 柄 


向 量 的 生成 


1) 直接 输入 向 量 
格式 上 要 求 向 量 元 素 需要 用 咱 " 括 起 来 ， 元 素 之 间 可 以 用 空格 、 过 号 或 分 号 分 隔 。 
用 空格 和 喜 号 分 隔 生 成 行 向 量 ， 用 分 号 生成 列 向 量 . 


2) 利用 冒号 生成 向 量 
冒号 表达 式 的 基本 形式 为 : X= Xx0:step:xn 
若 step = 1， 则 此 项 输入 可 以 忽略 。 


3) linspace 函 数 

可 以 使 用 linspace 函 数 生成 线性 等 分 向 量 : 

y=linspace(x1,Xx2) 生成 〈1 人 100 ) 维 行 向 量 ，y(1)=x1，y(100)=Xx2 
y=linspace(x1,x2, n) 生成 (1*n) 维 行 向 量 ，y(1)=x1，y(m=x2 


4) logspace 函 数 

logspace 用 于 生成 对 数 等 分 向 量 ， 格 式 如 下 : 

y=|logspace(x1,x2,n) 生成 (1*n ) 维 对 数 等 分 向 量 ，y(1)=10^x1， 
y(n)=10^Xx2; n 可 以 省 略 ， 此 时 其 默认 值 为 50。 








向 量 的 运 自 








人) 向量 加 减 与 数 加 减 


向 量 的 加 减 与 数 加 减 的 形式 与 普通 标量 加 减 相 同 
陪 人 向量 的 点 和 又 积 AS 一 混合 只 的 ) 实现 


点 积 : 向 量 的 点 积 函数 dot 实 现 。dot(a,b) 返 回 向 量 a 和 b 的 数量 点 积 ， 其 
中 a，b 兴 须 同 维 。 

又 积 : 又 积 由 cross 函数 实现 。 向 量 ab 必须 为 三 维 向 量 

混合 积 : 可 由 以 下 命令 实现 ，dot(a,cross(b,c)) 


3) 向 量 的 数 乘 、 数 组 乘 和 向 量 乘 
例 : 当 a= [1:1:3]; b=[2:2:6] 时 ， 以 下 命令 的 运行 结果 是 什么 ? 
2 2) >> a2=a.…b 3) >>ao=a*b 


,字符 串 : 包含 在 一 对 单 引号 中 的 字 





>> Ss='hello, MATLAB' % 定 义 字 符 串 变量 S 
S 一 
hello, MATLAB 
>> Whos 
NE ME 肖 的 ES 
S 1x13 26 char array 


Grand total is 13 elements using 26 bytes 


人 符 集 





单元 数组 〈(Cell Array ) CC 


单元 数组 是 MATLAB 数 组 的 一 种 特殊 数据 类 型 ， 它 用 于 保存 
不 同类 型 和 /或 不 同 大 小 的 数据 。 

三 种 直接 赋值 方式 

1， 单元 下 标 用 括号 " ( ) " 括 起 来 ， 而 单元 的 内 容 用 “ff" 括 起 来 ， 如 : 


>>Clear all 





2. 单 元 下 标 用 “" 括 起 来 ， 而 赋值 语句 等 式 右边 的 单元 内 容 用 中" 括 起 来 : 
>>alf1,1}=[1 2;3 4]; 

>>alf1,2}=[0 1 j; 

>>af2,1}='Hello';  % 右 边 只 有 一 个 元 素 时 可 省 略 去 "中 

>>al2,2}=2+93| 


3. 直接 使 用 人 
>>a={[1 2;3 4],[0 1], Hello,2+3i} 


单元 数组 的 操作 








>>a % 显示 单元 数组 a 的 信息 


先 使 用 也 数 cell() 创 建 空 的 单元 数组 ， 然 后 再 赋值 : 
>>b=cell(2,3) 

赋值 方法 同 直接 赋值 方式 。 

对 单元 数组 元 素 的 操作 

>>C=af1,2}】  % 将 单元 数组 a 的 {1,2} 元 素 赋 给 变量 C， 注 意 是 “ff ”而 不 是 


结构 体 区 沁 


。 与 C 语 言 类 似 ，MATLAB 结 构 体 用 于 存 取 相 关 的 
数据 ， 它 由 一 组 称 为 域 (fields ) 的 成 员 变 量 
(向 量 ) 构成 ， 每 一 个 域 可 以 为 不 同 的 MATLAB 
数据 类 型 。 
“结构 数 组 的 定义 有 两 种 方法 ， 一 种 是 直接 赋值 ， 
另 一 种 是 使 用 Strct() 邓 数 。 





结构 体 的 赋值 








>>Sstudent.name= Zhang Jun :; 

>>Sstudent.major=: Chemical Engineering ; 
>>Student.subject=[' 英 语 政治 ,数学 化工 原理 物理 化 学 ]; 
de] 性 =1811K21[e <》 《1 [7 ci 人 7 0 


>>Sstudent(2).name=Li Xia ; 

>>Sstudent(2).major='Chemical Engineering 
>>Student(2).Subject=[ 英 语 ”政治 ,数学 化工 原理 ”物理 化 学 站; 
>>Sstudent(2).entrance_exam=[60 72 68 85 88|]; 


struct_ array_name=structure(field1 ,values1,field2 ,values2,….) 
Student=struct(Cname ,Zhang Jun ，major, Chemical Engineering ) 





管道 压 降 的 计算 Ce 


T=[283:10:353]; 

M=2000;D=0.5; 
density.A=0.3471;density.B=0.274;density.Tc=647.13;denslty.n=0.28571 ; 
Rho=(density.A."density.B.^(-(1-T./density.Tc).Adensity.n))/0.2323; 
mu.A=-10.2158;:mu.B=1.7925e3:mu.C=1.773e-2;:mu.D=-1.2631e-5; 
mu=10.^A(mu.A+mu.B./T+mu.C<T+mu.D.T.^2); 

deltP=(M^1.8) (mu.^A0.2)./(20000"D^A4.8.*Rho) 


人 2 


也 数 文件 和 Script 文件 蛇 








Script 文件 
Script 仅仅 是 一 连 串 可 执行 的 MATLAB 命 令 ， 它 具有 全 局 性 
Script 文件 中 不 能 定义 函数 
函数 定义 的 一 般 格 式 : 
function [y1,y2,...,yn] = FuncName(x1;,Xx2,...,Xxn) % 函 数 声明 语 甸 


y1 一 ... 2%o (表达 式 1 ) 
y2 = ... % (表达 式 2 ) 
yn = ... o6 (表达 式 n ) 


其 中 ， 输 入 参数 为 X1,X2,...,xn, 输 出 参数 为 y1,y2,...,yn。 各 参数 可 以 是 
标量 、 向 量 或 矩阵 。 








邮 Editor -~ UntItJLed 因 








ile 了 dit TIext Cell Tools Debug Desktop 灿 indow Help 3 | 玉 居 
口 区 则 竹 品 芝 抽 号 下 交 和 目 怒 日 牛 Stack: 田 
四 四 | 


Script 于 [CLII 








阴 数 文件 的 编写 





编写 一 个 函数 ， 计 算 本 章 开 始 问题 中 流体 的 粘度 ， 
马 数 要 求 输出 在 度 的 计算 值 
. 邓 数 声明 语句 : function vis=viscosity() 
妆 量 的 传 交 
1. 通过 调用 函数 传递 
。 VIS=VviScosity(A,B,C,D,T) 
2. 通过 全 局 变量 传递 
。 利用 global 命 令 ， 在 主 函 数 和 子 函 数 中 予以 声明 
3. 编写 表达 式 


放 ”一 


画 数 的 调用 晤 





1. 在 调用 也 数 的 主 函 数 中 ， 直 接 采 用 函数 名 调用 
2. 通 过 马 数 铝 柄 调用 


YY 


管道 压 降 的 计算 函数 








eeliEGiiWEers ie 

eeleliLT 

=[25、 1 553| 

M=2000;D=0.5; 

density.A=0.3471; 

density.B=0.2740; 

density.Tc=647.13; 

density.n=0.28571 
Rho=(density.A.*density.B.^(-(1-T./density.Tc).^Adensity.m))0.2323 
mu.A=-10.2158;:mu.B=1.7925e3;:mu.C=1.773e-2;:mu.D=-1.2631e-5;: 
WESweSiN 

deltP=(M^1.8)”(Mu.^0.2)./(20000*DA^4.8.*Rho) 

5 人 
eteiNEEWEIeerSilN 

eeleliLT 

vis=10.^A(mu.A+mu.B./T+mu.C.“T+mu.D.“T.A2); 





朋 数 句 柄 Ge 


“创建 一 个 函数 多 柄 ， 可 用 于 保存 了 数 的 所 有 信 
息 ， 以 便 将 来 对 它 进 行 调 用 。 

。 马 数 名 柄 可 以 作为 站 数 ， 或 与 
feval 马 数 一 起 使 用 ， 以 调用 该 函数 铝 柄 所 属 的 函 
数 。 

。 使 用 函 数 包 柄 还 可 以 减少 定义 函数 的 文件 个 数 ， 
改善 重复 操作 的 性 能 ， 保 证 函数 计算 的 可 靠 性 

*funhandle = (Ofunction_name 
%function_name 为 用 尸 指 定 的 函数 名 


大 大 
1 
己 


道 压 降 的 计算 函数 





LeielE@iREeisiie2 

T=[283:10:353]; 

M=2000;D=0.5; 

density.A=0.3471 

density.B=0.274 

density.Tc=647.13; 

density.n=0.28571; 
Rho=(density.A."density.B.^(-(1-T./density.Tc).Adensity.n))/0.2323 
mu.A=-10.2158;:mu.B=1.7925e3:mu.C=1.773e-2;:mu.D=-1.2631e-5: 
WINMESweSiIVALLT NBGIE YESIweS IT 
deltP=(M^1.8)”(Mu.^A0.2)./(20000"DA^A4.8.*Rho) 


IleNESEY Eee 有 
vis=10.^(mu.A+mu.B./T+mu.C<T+mu.DT.^2); 





内 联 函 数 (inline function ) CC 
“ 内 联 函 数 是 Matlab 提 供 的 一 个 对 稍 ， 它 的 表现 和 函数 文件 
一 样 ， 但 内 联 函 数 的 创建 比较 容易 
。 内 联 函 数 的 创建 
vinline('CE)) 
vinline(CCE ,arg1,arg2,…) 
vinline(CE ,mn) 
。 涉 及 内 联 函 数 性 质 的 指令 
wclass(inline_fun) 同和 
wcharkinline_fum) 给 出 内 联 函 数 计 算 公 式 
vargnames(inline_fun) 给 出 内 联 函 数 的 输入 变量 
vvectorizelinline_fun) 使 内 联 函 数 适 用 于 数组 运算 
。Matlab 中 许多 "“ 泛 了 ? 画 数 都 是 采用 inline， 从 而 具备 了 适应 
各 种 被 处 理 函 数 形 式 的 能 








内 联 双 效 的 应 用 


动因 [IIS 人 《ITIe7aaie 
f{=F1(2) 


站 二 etelipxs( 
XX=[0.5,1,1.5,2];f1=FF1T(XX) 


G2=inline('a*exp(x(1))*cos(X(2))，a,X) 
g1=G2(2,[-1,pi/3]) 





匡 芝 人 EeeSainiiwie CC 


。 匿 名 函数 用 于 在 命令 行 、 函 数 文件 或 script 文 件 中 创建 简 
单 形式 的 函数 ， 避 免 另 外 定义 新 的 函 数 
。 匿 名 马 数 的 定义 形式 
避 CItillelSI 检 esSiel 


f@(x) x.^2 人 人 @(X) X.^2: alpha=0.9: 
a=f(5) g=@(X) 3 Xi f-@(x) sin(alpha*x); 
结果 : a= 25 gf f(pi) 


结果 : ans = 0.3090 


数据 输入 和 输出 





。 数 据 输入 
一 利用 M 文 件 产 生 数据 文件 
-用 Load 人 驹 令 从 MAT 文 件 或 文本 文件 读 取 数据 
一 用 fScanf 了 数 
-用 提示 输入 函数 input 
二 else 关 [jeelite[tc 本 人 SteiEE 
。 数 据 输出 
-用 Save 命 令 
一 用 fprintf 驰 数 
-用 函数 disp() 将 结果 输出 至 屏幕 
-dlImwirte,xlswrite 马 数 
一 图 形 输出 


Matlab 二 维 图 形 


数据 准备 : 

e 选 定 所 要 表现 的 范围 

e 产 生 自 变量 采样 向 量 

@ 计 算 相 应 的 函数 值 回 量 

选 定 图 形 窗 及 子 图 位 置 

@ 缺 当时， 打开 Figure No.1， 或 当 六 窗 ， 
当前 子 图 

@ 可 用 指令 指定 图 形 窗 号 和 子 图 号 
调用 《局 层 ) 绘图 指令 ; 在 指令 中 设置 
线 型 、 色 彩 、 数 据点 型 

设置 轴 的 范围 与 刻度 、 坐 标 分 格 线 














图 形 注 酝 : 
图 名 、 坐 标 名 、 图 例 、 文 学 说 明 





图 形 的 精细 修饰 《图 柄 操作 ) : 
@ 利 用 对 象 属性 值 进行 设置 
@ 利 用 图 形 窗 工 具 条 进行 


Ce 


t=pix(0:100)/100; 
y=SIin(t).xSin(9K#t; 


figure(D) 上 %% 指 定 1 号 网 形 窗 
subplot(2,2,3) %% 指 定 一 个 具有 2 行 2 列子 
图 形 窗 中 的 3 号 子 图 


plot(ty,b-) “多 用 蓝 色 实 线 画 曲线 


axis([0,pi,-1,1]) 多 设置 轴 的 范围 
grid on %% 男 举 标 分 格 线 


title(' 调 制 疲 形 ) 多 图 名 
xlabel(Ct);ylabel(y) 双 轴 名 
TANGITONNIVONITG 六 帮 72 
text(2,0.5,y=sin(bsin(9D) 狗 文 字 说 明 


set(h,'MarkerSize',10) 双 设 置 数据 点 大 小 





函 数 Plot 基 本 调用 格式 3 








1) plot(X,'S)) 

。 X 为 实 向 量 时 ， 以 该 向 量 元 素 的 下 标 为 横 坐 标 、 元 素 值 为 纵 坐 标 画 一 条 连续 曲 
线 

。 X 是 实 和 矩阵 时 ， 则 按 列 绘制 每 列 元 素 值 相对 其 下 标的 曲线 ， 曲 线 数 目 等 于 X 的 
列 数 


。 X 是 复数 天 阵 时 ， 则 按 列 分 别 以 元 素 的 实 部 和 上 庶 部 为 横 、 纵 坐标 绘制 多 条 曲线 
。 'S' 是 用 来 控制 线 型 、 色 彩 、 数 据点 型 的 选项 字符 串 。S 可 以 缺 省 ， 此 时 曲线 按 
Matlab 默 认 设 置 绘制 。S 的 取 值 见 下 节 

2) plot(X,Y,'S) 
。 X、Y 是 同 维 向 量 时 ， 绘 制 以 X、Y 为 横 、 纵 坐标 的 曲线 
。 X 是 向 量 ，Y 是 有 一 维 与 X 同 维 的 和 矩阵 时 ， 则 绘 出 多 根 不 同色 彩 的 曲线 。 曲 线 
数 等 于 Y 的 另 一 维 ，X 作 为 这 些 曲线 共同 的 横 坐 标 
。 X 是 和 矩阵，Y 是 向 量 时 ， 情 况 与 上 相同 ， 只 是 曲线 都 以 Y 为 共同 纵 坐 标 
。 X、Y 是 同 维 矩阵 时 ， 则 以 X、Y 对 应 列 元 素 为 横 、 纵 坐标 分 别 绘制 曲线 ， 曲 线 
条 数 等 于 和 珑 阵 的 列 数 。 
。 'S' 的 意义 ， 与 上 相同 。 

3) plot(X1,Y1,'S1',X2,Y2,'S2.….) 
。 此 格式 中 ， 每 个 绘 线 “ 三 元 组 ”( X,Y,'S' ) 的 结构 和 作用 ， 与 上 相同 。 不 同 “ 三 元 
组 ”之 间 没 有 约束 关系 。 


曲线 的 色彩 、 线 型 和 数据 点 型 貌 





条 号 
_- 和 包 吕 醒 避 本 
线 


AAA 符号 | 


符号 | 


IE 








坐标 、 刻 度 和 分 格 线 控制 Ge 





本 T 到 








axis auto 缺 省 设置 纵横 轴 采 用 等 长 刻度 
坐标 充满 整个 绘图 区 
axis off 取消 轴 背 景 ET 纵 、 横 轴 采 用 等 长 刻度 ， 
| 坐标 框 紧 贴 数 据 范围 
TIETTTT 本 
axis j 矩阵 式 坐 标 ， 原 点 在 左上 产生 正方 形 坐 标 系 
方 
axis xy 普通 直角 坐标 把 数据 范围 设 为 坐标 范围 
axis(V) 人 工 设 定 坐 标 范围 。 设 定 | 。 axis vis3d 保持 高 宽 比 不 变 
WSIR2w5 2 站 有 二 省 外 及 
三 | 芭 上 xx 攻 V2 





诈 导 M 


图 形 标 识 


图 形 标 识 命令 调用 格 却 


EU 二 TITLE(text) 
LiejeiN eei=NALLLIL ee 
ropertyValue2,…) 





从 标 机 (label ) XLABEL(Ctext ) 
XLABEL(text,Propertyl1 ,PropertyValue1,Property2 ,Pr 


NT 


和 图例 〈legend ) LEGENDI(string1,string2,string3，.…) 








名 是 ce 


T=[283:10:353]; 

M=2000;D=0.5; 
density.A=0.3471;density.B=0.274;density.Tc=647.13;density.n=0.28571; 
Rho=(density.A."density.B.^(-(1-T./density.Tc).^Adensity.n))/0.2323; 
mu.A=-10.2158;:mu.B=1.7925e3;mu.C=1.773e-2:mu.D=-1.2631e-5: 
mu=10.^A(mu.A+mu.B./T+mu.CT+mu.D.T.^2); 

deltP=(M^1.8) mu.^0.2)./(20000-D^4.8…Rhoj)plot(T ,deltP,b-o)) 

title( The pressure drop vs Temperature ') 

xlabel(Temperature(^oC)) 

dlelsl 人 adsisisigliselteje 芭 (eelllN《ziL[=]81 只 [= 下 e1 首 ejjeis 翅 


陋 必 


古国 和 
| 和 


和 和 fontnatmefatrgl | 且 有 windowrs | %fontnatmefcourieryEzatmnhlel， | Examplel 
字库 字体 | fontname( 隶 书 ) 范 网 2 范 团 2 


负 记 Batmnhle 了 EXamp|e 全 
ASatT1Ed4 后 Y 呈 后 折 


和 要 十 EL 二 } 下 Batr 站 1 
Di 二 ET 站 了 蕊 Batm 和 站 


Example5 


EnPle 癌 





句柄 图 形 的 结构 CS 





句柄 图 形 的 访问 与 操作 




















依 属 性 值得 找 铬 形 对 象 
获得 对 象 属性 


设置 对 象 属性 





删除 儿 形 对 象 
几 信 








例题 


绘制 二 阶 系统 阶 跃 响 应 曲线 





例题 区 








clf;t=6*pi*(0:100)/100;y=1-exp(-0.3*t).*cos(0.7*t); 
tt=t(find(abs(y-1)>0.05)); 9% 了 寻找 大 于 0.05 的 元 素 
ts=max(tt); % 寻 找 tt 中 最 大 的 元 素 
plot(t,y,r- LineWidth ,3) 

axis([-inf,6*pi,0.6,in 人 ]) 

set(gca, Xtick',[2*pi;,4*pi,6*pi],，Ytick'",[0.95,1,1.05,max(y)]) 

grid on 

title(Nity = 1 - e^A{f -alphatjcosf{yomegat} ) 

text(13.5,1.2, \fontsize{12}falpha}=0.3)1) 

text(13.5,1.1, \fontsize{12j 人 omegalj=0.7”) 

hold on;plot(ts,0.95, bo' ,MarkerSize ;10);hold off 

cell_string{1}= \fontsize{12 和 uparrow ; 

cell_string{2}=\fontsize{16} \fontname{ 隶 书 } 镇 定时 间 "， 
cell_string{3}=\fontsize{6} “; 

cell_string{4}=[\fontsize{14})\rmt_ {s} = " num2strlts)]; 
textlts,0.85,cell_string) 

xlabel(\fontsize{14} \bft \rightarrow ) 

EL 人 NIeliGir43 LN LEULeTY 人 





双 。 六 AT 
多 次 王 绘 、 双 纵 坐 标 和 多 子 图 


人 多 次 滞 绘 : 


hold on， 使 当前 轴 及 图 形 保持 而 不 被 刷新 ; hold off， 不 保持 当前 轴 及 图 形 。 


2) 双 坐 标 图 : 

plotyy(X1,Y1,X2,Y2) 以 左右 不 同 纵 轴 分 别 绘制 X1-Y1、X2-Y2 两 条 曲线 
plotyy(X1,Y1,X2,Y2,Fun) 以 左右 不 同 纵 轴 把 X1-Y1、X2-Y2 绘 制 成 Fun 指 定 的 形式 
的 两 条 曲线 


3) 多 子 图 : 

采用 subplot(m,n,k) 使 (m xn) 幅 子 图 中 的 第 k 个 成 为 当前 子 图 ， 再 采用 其 它 的 图 形 
绘制 指令 则 可 将 图 形 绘制 到 指定 的 子 图 中 。 子 图 序号 的 编制 原则 是 : 左上 方 为 第 1 
幅 ， 向 右 向 下 依次 增 大 。 





双 坐 标 曲 线 绘制 方法 cG9 
吾 出 函数 ET 国生 入 分 国人 和 刘 在 区 辣 加 四 上 的 
曲线 ， 





clf; 
dx=0.1;x=0:dx:4;y=X.*sin(X);jS=cumtrapz(y)*dxX; % 梯 形 法 求 累 计 积分 
plotyy(X,y,X,S),text(0.5,0, \fontsize{14Nity=XxsinX ) 

Sint= 人 fontsize{16Nint_ ffontsize{8}0}A{ X]}) 
text(2.5,3.5,[\fontsize{14》Nits= ,sint, \fontsize{14NitxsinXxdx ]) 


Matlab 三 维 图 形 








1) 三 维 曲线 绘制 镍 令 plot3 
2) 三 维 网 格 图 形 绘制 命令 mesh 
3) 三 维 曲面 绘制 名 利 sur 


[X0O Y0 Z0]=sphere(30); 
X=2*X0;Y=2*Y0;Z=2*2Z0; 
Surf(X0,Y0,Z0) 
tellileailtcie 

Jelegell 

mesh(X,Y,Z) 
welielduiiertiteiy 
[ee 

[elleEelii 

axis equal 

axis off 





Matlab 图 形 绘 制 函 数 


Matlab 图 形 绘制 函数 分 属于 以 下 帮助 主题 
。 graphe2d 
ellellkeie 
je 过 weltiei 





程序 的 组 成 





Script 文件 或 国 数 文件 
RCR 
数据 输入 3 运 入 一 数据 竹 出 
键 嫩 输入 数学 运算 碾 腊 竹 出 
文件 输入 关系 运算 文件 竹 出 
也 和 辑 运算 图 形 得 出 


流程 控制 


MATLAB 关 系 运 算 符 











MATLAB 逻 辑 、 关 系 运 算 的 规定 CS 


。 交 辑 和 关系 运算 均 可 以 在 任意 具有 相同 维 数 的 数组 之 
间 进 行 

。 标量 可 以 和 任意 维 数 的 数组 运算 

。 在 所 有 关系 表达 式 和 逻辑 表达 式 中 ， 作 为 输入 的 任何 
非 0 数 均 被 看 作 是 “逻辑 真 "”， 只 有 0 是 “逻辑 假 ” 

。 所 有 关系 表达 式 和 还 辑 表达 式 的 运算 结果 是 一 个 由 0 
和 1 组 成 的 膛 辑 数组 

。 浆 辑 数组 是 “数值 数组 ?的 子 类 ， 可 以 按 数值 数组 实施 
操作 ; 它 又 有 不 同 与 普通 数值 数组 的 特性 ， 即 表示 着 
对 对 象 判 断 的 真 假 ; 可 用 于 数组 寻访 等 特殊 场合 


MATLAB 巡 上 辑 运 算 符 








A-=[0110]; B-[1100] 


Ag&B-[0100] AIB-I1110] ~A=I1001] xor(AB)-=I1010] 
先决 还 辑 运 算 符 : 


先 关 或 | 1 


次 辑 、 关 系 运算 示例 


A=-3:2; 
B=3:-1:-3; 
L1=~(A>0) 
L2=~A>0 
L3=~A 
L4=A>-2&A<1 
L5=A==B&L1 


吕 口 口 口 一 


加 一 口 口 一 


本 


CD GD 所 


CE De 一 


口 口 口 口 口 上 品 








运算 优先 级 : 


括号 癌 
逻辑 否 呆 
乘除 
加 减 

关系 运算 
逻辑 运算 
先决 与 
先决 否 


MATLAB 关 系 运 算 函 数 


ISempty 
己 >elgl 
any 

all 

le 
SSeczilU 
[Sewrteli 
[1 
ISInf 
Siillte 


数组 是 否 为 空 

两 个 数组 是 否 相 等 

数组 有 非 零 元 素 则 结果 为 1 
数组 元 素 全 非 零 则 结果 为 1 
数组 非 零 元 素 的 下 标 
是 否 为 标量 

是 否 为 向 量 

是 否 为 非 数 

是 否 为 无 穷 


是 否 为 有 限 








for 循 环 结 构 


for 衢 环 变量 = 表达 式 1 ( 初 值 ) : 表达 式 2 ( 步 长 ) : 表达 式 3 ( 终 值 ) 
statements ( 语 负 组) 
end 字 符 串 : 包含 在 一 对 单 引 号 中 的 字符 集合 


for i=1:10 
X(I)=|; 
ie 

X 


。 为 了 得 到 局 效 代码 ， 应 尽量 提高 代码 的 癌 量 化 程度 ， 避 免 使 用 循环 
结构 


。 为 了 得 到 局 效 代码 ， 在 循环 指令 之 前 应 尽量 对 数组 进行 预定 义 





while 循 环 结构 


while cndition (表达 式 ) 
statements ( 执行 语 名 组 ) 
ie 


Fibonacci 数 组 的 元 素 满足 Fibonacci 规 则 : 
alkiz=ak+ak:1，〈k=1,2,…); 且 al1=a2=1。 求 该 数组 中 第 一 
个 大 于 10000 的 元 率 。 





if-else-end 分 支 结构 








if 语 铝 的 一 般 格式 : 
| 芭 welllelillelill 
statements 。 % 如 果 condition1 的 值 为 True， 则 执行 该 语 多 组 
= 已] 攻 welilellil[eli 挛 
statements 。 % 如 果 condition2 的 值 为 True， 则 执行 该 语 负 组 
else 
statements 。 % 如 果 condition1 和 condition2 的 都 为 False， 则 执行 该 语句 组 
ie 


用 for 循 环 指令 寻求 Fibonacci 数 组 中 第 一 个 大 于 10000 的 元 素 
n=100;a=1:100; 
for I=3:n 
a(l)=ali-1)+al(i-2); 
f a()>=10000 
al 人 由 ， 
break; 
ie 
ie 
| 








Switch-case 结 构 








Switch-case 的 一 般 格 式 : 
Switch test_expr % 测 试 表 达 式 test_expr 可 以 是 标量 或 字符 串 
case value 
statements ”%% 当 test_expr 值 是 value 时 ， 执 行 该 语 多 组 
case {value1,value2，….} 
statements “”% 当 test_expr 值 是 value1 或 value2 或 .………… 时 ， 执 行 该 语 多 组 
otherwise， 
titsTsalS 
ie 


Switch...case...otherwise 语 可 的 能 力 与 并 ...else...end 语 名 类 似 ， 
但 对 多 重 选择 的 情况 switch 语 名 使 代码 更 加 易 读 。 


try-catch 结 构 








try-catch 的 一 般 格 式 : 
try 


statements “”%% 此 组 语句 总 被 执行 。 知 正确 则 味 出 此 结构 
catch 


statements “”% 当 上 组 语 名 出 现 执行 错 误 后 ， 该 组 语句 被 执行 


门 
站 N-=4;A=[1 2 3] 
try 
A_N=A(N) 

。 当 两 组 语句 都 出 错 后 ，Matlab 将 catch 

跳出 该 结构 A_N=Alend) 

。 可 以 采用 lasterr 马 数 查询 出 错 原 end 
lasterr 
AN- 3 
amns = 三 


Index exceeds matrlx dimensions. 





本 讲 小 结 


RE 算 有 函 数 的 使 用 
。 注意 区 别 和 天 阵 和 数组 的 乘 、 除 、 乘 方 运算 
atlab 数 据 输入 输出 功能 ， 克基 是 从 会 图 功能 的 实现 
atlab 马 数 文 件 的 基本 形式 及 其 调用 
字符 、 单 元 数组 、 结 构 体 的 定义 
5)Matlab 的 流程 控制 语 多 


第 二 讲 
非 线性 方程 (组 ) 求解 与 迭代 法 


化 工学 院 软件 应 用 教科 组 
2006-10 


本 章 知 识 要 点 





。 数值 计算 
- 单个 非 线性 方程 求解 
- 非 线性 方程 组 求解 
- 迭代 法 


。 MATLAB 
- 求解 非 线性 方程 《组 ) 的 相关 函数 


本 章 的 所 要 解决 的 典型 问题 
在 945.36kPa (9.33atm ) 、300.2K 时 ， 容 器 中 充 以 
2mol 气 气 ， 试 求 容器 体积 。 
已 知 此 状态 下 氮气 的 P-V-T 关 系 符合 范 德 华 方程 ， 其 
范 德 华 常 数 为 a=4.17atm*L/mol2，b = 0.0371L/mol 








多 组 分 混合 浴 液 的 沸点 、 饱 和 蒸气 压 计算 
* 流体 在 管道 中 阻力 计 租 

“多 组 分 多 平衡 级 分 离 操作 模拟 计算 
| 数 法 求解 化 学 平衡 问题 

“ 定 态 操作 的 全 混流 反应 器 的 操作 分 析 


求解 非 线性 方程 的 方法 





。 二 分 法 

。 不 动 点 迭代 

。 威 格 斯 坦 法 欠 代 
。 午 顿 法 

。 钊 | 线 法 





认 代 法 原理 


0 = 0 一 = (00 
己 知 : -xz-1=0| 


的 解 为 1.324， 以 1.5 
为 初 值 ， 采 用 以 下 
网 种 迭代 格式 计 

算 ， 结 来 如 何 ? 





和 迭代 法 意义 示意 图 








X1 X3 X2 
-1< 中 "00<0 


X2 X1 X3 
中 "OO<-1 





X3 X2 XI1 





0) = 
Matlab 表 达 多 项 式 方法 ， ET 和 





Matlab 多 项 式 了 数 
- Polyval 多 项 式 的 值 
- Polyder 多 项 式微 分 
- Polyfit 多 项 式 拟 合 
ie 人 BDI=ieelii 多 项 式 乘 法 和 除法 


-Roots 多 项 式 来 根 


简单 多 项 式 函 数 的 使 用 


a=[1 11 55 125]; 
b=[1 1;1 1j; 
eelelNA《L[( 有 ie 
q=[2,3,5]: 

te 世 well(Rej 

e 园 elelNAALGLej 
= 二 eeeliN(cUeE2j 





% 求 多 项 式 a 在 b 处 的 值 


% 多项式 乘法 
% 多 项 式 向 量 表 示 为 符号 多 项 式 
% 多 项 式 除 法 


多 项 式 求 根 函 数 roots 
求解 方程 : 


ie 司 Leleitife) 
结 采 : 


SO| = 


1.6024 + 1.27091 
1.6024 - 1.2709 
-0.3524 + 0.9759| 
-0.3524 - 0.975951 








PE 9 .3.3， 

= 300.2; 

n = 2: 

a = 4.17: 

b = 0.0371; 

R = 0.08206 

coef=[P, -(P nb+mR TD an^2 ,- 
ab-n^3]; 

人 elIelSi(eels il 

结 采 : 

V =5.0028 


0.24<9 “0.1092 





非 线性 方程 来 解 函 数 fzero 


调用 格式 : 
[x,fval,exltflag,outputl = fzerolfun,x0,options, p1，pP2，.….) 
此 函数 的 作用 人 的 零 值 点 ，X0 是 标量 。 


《zl 数 在 解 X 处 的 值 

eXxitflag 3 直 束 情况 ，>0， 程 序 收 敛 于 解 ; 
<0， 程 序 没 有 收敛 ; =0， 计 算 达 到 
了 最 大 次 数 

one 一 个 结构 体 ， 提 供 程 序 运行 的 信息 ; 


output,iterations， 迁 代 次 数 ; 
output.functions， 邓 数 fun 的 计算 次 数 ; 
output.algorithm， 使 用 的 算法 
ejeliieli 选项 ， 可 用 iopimeaiz 数 设 定 选 项 的 新 值 
fun 可 以 是 函数 句柄 避 匿 名 函数 


fzero 亏 数 的 使 用 


1) sinXx 在 3 附近 的 零点 


2) coSsXx 在 [1,2] 范 围 内 的 零点 


-2sinz=0 


1 


PileI(CjIRG) 
fzero(@2cos,[1,2]) 


fzero(@(X) X^A3-2*Ssin(X),1) 
或 者 : 
LewieEeiit 志 eslilel 
X=fzero(fun,1) 
ee 

y=X^A3-2 -Sin(X); 


fzero(@(x) X^A3-2*X-5, 人 ); 
roots([1 0 -2 -5]) 








fzero 马 数 初 值 的 选取 





民 二 日 J3， 

= 300.2; 

n = 2: 

a=4.17; 

b = 0.0371; 

R = 0.08206; 

V=fzero(@(V) PVA^3- 
(Pnb+rR TD 六 VA2+ an^A2V - 





fzero 函 数 初 值 的 选取 
7JO=Gn De 一 川 县 于 汪 


以 t 为 自 变量 ， 取 值 范 围 为 -10<t<K10，ab 为 参数 ， 本 例 取 值 分 
别 为 0.1，0.5 


LeleliE@iiLT 二 elsiilie 产 
a=0.1;b=0.5;t=-10:0.01:10; 
Y=sin(t).^2.*exp(-a*t)-b*abs 人 tb); 
9 前 el et 的 的 的 用 elie 区 el1i 机 elieitf 罗 XierSi(Sjpxk0) 风 全 
[1e1=] [外 时 区 1[21e=] 的 《 人 二 册 eliegelil 
zoom on 
n=input(CHow many zero points are there?) 
本 4 蔬 et)pxelelilEelil 
for 1=13:n 

[to(),y(),exitflag]j=fzero(@ 由 sin(t)^2*exp(-a*t)-b*abs(t),tt(iD)); 
ie 
disp( The zero points are2) 
fprintf(%e.4ft ,tO) 
fprintf( An ”) 


fzero 马 数 初 值 的 选取 








在 945.36kPa (9.33atm ) 、300.2K 时 ， 容 器 中 充 以 2mol 气 
“本 试 求 容 器 体积 入 。 

已 知 此 状态 下 气 气 的 P-V-T 关 系 符 合 范 德 华 方 程 ， 其 范 德 华 时 
0 L/mol2，b = 0.037 并 /0 


function PVT 

P = 9.33; ?%e atm 

T = 300.2; % K 

n=2; 9%omol 

a = 4.17: 

b = 0.0371 ; 

R = 0.08206; 

V0 = RTV/P 

品 fval] = fzero(@PVTeqV0, 由 ,P,T,n,a,b,R) 


function f= PVTedq(V,P,T,n,a,b,R) 
f=(P+asn^2VA^A2) ”(V-nb) - mmR Ti 


预 热 到 了 的 会 有 尽 放 物 的 谊 祝 原 料 ， 以 一 定 的 六 量 只， 加 六 到 容 可 汶 三 的 搅拌 
模 民 旗 研 中 进行 绝热 民 应 。 尽 应 福全 物 连 续 排 出 。 衬 的 进 、 出 口 六 度 分 别 济 Can 
和 上。 尽 放 带 袍 的 密度 光 PP， 比 热 雁 为 cp。 模 内 及 出 口 刘 度 汶 T， 民 旗 速 度 汶 ， 
-Pr=C3， 忒 中 直 = 丰 egb 全 二) 。 已 郑 数 据 ， To=450 及 ，Cm 人 吕 )AoC， = 250K， 


ER=10000R， 大 它 em， 工 = 防 / 让 = 总 2， 试 求 反应 转化 率 。 
避 型 ， 由 物料 往 算 和 枚 量 简 章 可 以 获 行 蛋 型 方 杜 如 下 
bcatl- aremp 人 -xz=0 (物料 衡 算式 ) 


7- 了 = 王 2 5 【热量 衡 算式 ) 


电 


LewlelEGeN 人 iT 汐 elsiiile7 
T0 = 450; 
x0 = [0:0.12:1.0]; 
n=|length(XxO); 
for 1I=1:n 

x()= fzero(@NonlinEq,x0() 咯 ,TO); 
iie 
[je 全 由 Leellh《sltsileliEewelillleale 
fprintf(%e.4fNt ,XI) 
fprintf( An ) 
人 
Leile]iiE Ne 过 er 和 全 人 0) 
T=T0O+250”X; 
f{ = 0.25"(1-X)^A2-exp(20-10000/T) - X; 


结 采 : 

esISi[eIiEweill[egleIs) 

0.0408 0.0408 0.2804 0.2804 0.2804 0.2804 
0.8363 0.8363 0.8363 0.8363 1.0951 





fsolve 函 数 


。 与 fzero 函 数 只 能 求解 单个 方程 的 根 不 同 ，fsolve 邓 数 可 求 
解 非 线性 方程 组 的 解 。 其 算法 采用 的 是 最 小 二 乘法 。 
。 调用 格式 : 

| 人 NI 区 =》>41i[Ue 区 elgliiel 员 晤 [=Uweleltc16] 攻 二 二 [ielN《sI4ITL 了 《6 下 eljeltlelaT 汉 
p1, p2,，.…) 
。 输入 输入 变量 的 意义 同 fzero 函 数 。 输 出 变量 中 的 jacobian 
为 函数 fun 在 x 处 的 Jacobian 拢 阵 。 


fsolve 马 数 的 使 用 








LatelEG@iLT2eisiiliels 
Xx0O=[1 1 1j; 
X=fsolve((Ofun,X0O) 
function y=fun(X) 
y(1)=SIn(X(1))+X(2)^A2+Iog(X(3 放 -7 
y(2)=3*X(1)+2^Xx(2)-X(3)^A3+1; 
Yy(3)=X(1)+X(2)+X(3)-5; 

解 得 结果 如 下 : x 王 0.5991 ”2.3959 2.0050 





fsolve 马 数 的 应 用 


在 铜 管内 在 1 atm 下 将 异 丙 醇 加 热 到 400"C 。 已 知 铜 是 生产 丙酮 和 丙 醛 的 
催化 剂 ， 或 许 还 有 某 些 异 两 醇 异 构 化 为 正 丙 醇 。 这 三 种 产物 的 生成 可 用 
如 下 三 个 独立 反应 表示 : 

ex] re 人 (| 二 用 exlarela[ 人 | 2 K1 = 0.064 

ee (| 二 三 福 (e kjIeelf GT: 

iC3H7OH(IP) =C2H5 CHO(PR) +H2 K3 = 0.00012 

后 续 加 工 步 又 需要 正 丙 醇 ， 虽 然 可 含 丙 酮 ， 但 两 醛 含 量 不 能 超过 
0.05(mol)%e。 在 上 述 反 应 条 件 下 ， 是 否 存在 违反 这 种 规定 的 可 能 性 ? 
数学 模型 : 各 反应 的 化 学 平衡 方程 如 下 











LedeE@iii 记 ee 
= [0.05 0.2 0.01]; 
x = fsSolve(@EquiC3,x0); 
CAC=x(3)/Sum(X) 
f CAC<0.05 
[ie 全 中 LV@XeelilwsiilltztielEeelllieglitelmleise= 训 0 下 65 
> 
ee 全 有 [weliiteslitgzileiEeellieglei se 帮 0 史 0 
ie 
function f = EquiC3(X) 
f1 = Xx(1)-0.064*(1-X(1)-X(2)-X(3)); 
f2 = X(2)*(X(2)+X(3))-0.076*(1-X(1T)-X(2)-X(3))*(1+X(2)+X(3)); 
f3 = Xx(3)*(X(2)+X(3))-0.00012*(1-Xx(1T)-X(2)-X(3)) (1+X(2)+X(3)); 
f= [ff2f3]; 


welitesTiligztiieii 区 eeigl[eElirelmlei= 太 es 0 页 0 从 


水 平 管道 内 晖 定律 流体 的 闹 动 - 非 牛 顿 访 体 闹 动 所 需 管 径 

一 管道 内 咎 定 刍 流 体 的 水 平流 动 ， p= 964kgim3， 质 量 流量 GE667kg's1， 管 长 
T=10im，gE= 5X105m，ap= 15kPa 攻 =18，n= 妈 ,及 = 148M satm2。 试 用 
普遍 作 的 和 雷诺 数 法 估算 此 非 牛 顿 入 体 访 动 情况 所 需 的 管道 直径 。 

数学 模型 ， 





普 这 化 的 雷诺 用 法 定义 ，Re us = 富 -Q ， 式 中 ，a 和 尺 是 前 切 速率 的 函 台 


对 于 埋 定 律 流体 ，n= nm， 天 - 对 世 生 | 





2 
普遍 化 的 Re 二 系 也 : 
层 流 时 (Few <2100) 7 了 = 二 
人 
渗 访 时 【Re > 2100 厅 =-7te Ba Re > 


委 道 直径 D 的 计算 是 选 从 过 程 ， 其 起 冶 鸽 计 值 按 下 式 估 看， 


,区 38+11 
32 开 8 4CY 
2 人 | ”On 
闻 TB 


站 一 





RE 加 
D 的 改进 估计 值 ， 局 2 6 .其 中 心 = 各 


DLReiEPeLe 
rho = 961;G = 6.67jiL = 10;eps = 5e-6;deltaP = 15e3;k = 1.8;n = 0.64;K = 1.48;n1 = ni 
K1 = K*((3*n+1)/(4*n))Ani 
D1 = (32*K1*L*8^(n1-1) (4*G/pirho)An1)A(1/(3*n1+1)); 
D2 = 2*D1; 
delta = 1e-4; 
while abs((D2 - D1)/D2) > delta 
uU = (G/rho)/(pFD1^A2/4); 
Regen = D1^n1*uA(2-n1)*rho/(8^(n1-1)*K1); 
if Regen <= 2100 
f= 16/Regen; 
end 
if Regen > 2100 
f= FrictFactor(Regen,n1); 
| 
Le = k*D1/(4*f); 
D1 = D2; 
D2 = (2*f(L+Le)/(rho*deltaP)*(4*G/pi)^2)^0.2; 
end 
fprintf("t 管 道 直 径 为 D = %.4f %s\n',D2,"m'") 
fprintf("\t 摩 擦 因子 为 f = %.4f %s',fm') 
9 ----=-----=-======= -==== 一 = 一 = 一 === 一 = 一 = 一 == 一 一 = 一 一 一 
function f = FrictFactor(Regen,n1) 
f0 = 16/Regen;i ”9%4.53 
f = fzero(@fFunc,f0,[,Regen,n1) 
0 全 
function F = ffunc(f) 
F = 1/sqrt(f) - 4.0*log(Regen*f^(1-n1/2))/n1^0.75 + 0.4/n1^1.2;2%64.54 


本 讲 小 结 





Let 


w [sol,feval,exitflag,outputj=fzero(Cfun， 
X0;,options,p1,p2,.….) 


:AIR LEU eeI EUweeicIDEie 
MI(CIIIIRB LI ITe IE 因 oP24 交 


第 三 讲 
答 阵 操作 与 线性 方程 组 求解 


化 工学 院 软 件 应 用 教科 组 
2006-10 


本 章 知 识 要 点 





。 数 信 计 复 
- 线性 方程 组 的 直接 解法 
- 线性 方程 组 的 迭代 解法 





。 MATLAB 
- 宅 阵 的 生成 与 操作 
- 窍 阵 运算 函数 











7% 二 甲 薪 
49%6 葵 乙 燃 


549% 甲 茉 
359% 茉 
18% 二 甲 蒜 
159% 二 249% 苯 乙 烯 
25% 。， ] 429%6 甲苯 
409% 甲 葵 16% 苯 
20% 匠 D2 15% 二 甲 蒜 
4 10% 蒜 乙 烯 
F=7o kg-molAmin 5496 甲 某 
219%6 茉 


24% 二 甲 茉 
659% 葵 乙 烯 
10% 甲 茱 
19%6 匠 


























B2 


>》(FGn) (CD)= >》(F(CoxD (ou1) 


线性 方程 组 在 化 工 模型 中 的 作用 





。 多 组 分 体系 的 物料 衡 算 
。 各 种 化 合 物 的 物理 化 学 性 质 
。 稳 态 动力 学 计算 

。 微 分 方程 的 差分 方程 求解 


病态 线性 方程 组 


X=1,y=1 x=10,y=-2 


x=0.991,y=-0.487 








X=2,y=-2 


线性 方程 组 的 求解 方法 





十 接 求解 方法 : 
。 局 斯 消 元 法 
。 局 斯 一 约 当 消 去 法 
。 扔 赶 法 
人 适 代 解 法 : 
-- 雅 可 比 (Jacobian ) 迭代 
- 高 斯 - 赛 德尔 (Gauss- Siedel ) 迭代 
一 松弛 (SOR ) 迭代 
-- 共 力 梯度 (CG ) 和 迭代 
在 Matlab 线 性 方程 组 的 解法 归结 为 矩阵 的 除法 命令 


短 阵 的 生成 


直接 输入 小 矩阵 
了 人 二 elire] 
“创建 M 文 件 输入 大 短 阵 
“利用 天 阵 操作 命令 
一 Cat - 数组 连接 
-reshape -数组 变 维 
本 11 





常用 短 阵 生成 函数 





形 有 具 憩 此 ie 上 - 全 零 阵 
Js - 全 一 阵 
eye -单位 阵 
EU - 均匀 分 布 随机 阵 
Te - 正 态 分 布 随机 阵 
B=zeros(mn) rand(state ,mn) 
B=zeros(m,n) rand(N) 
B=zeros(Slize(A)) 
State 
twister 


Seed 


特殊 珑 阵 


常用 短 阵 生成 函数 


weliljei1l 
gallery 


teEUEUe 


hankel 
hilb 
Le 
Etelle 
ee 
[esi 
toeplitz 
We[] 
WMS 








ieliile1ilieliiliitzttih4 

- Higham test matrices. 
Eeegiiictthe 

- Hankel matriX. 

- Hilbert matriX. 

- Inverse Hilbert matriX. 
EUel[eisiell1t 

- Pascal matriX. 

- Classic symmetric eigenvalue test problem. 
ee 

-Vandermonde matriX. 

- Wilkinson's eligenvalue test matriX. 





一 维 数组 的 寻访 和 赋值 


生成 一 (1x5) 0 状态 下 的 均匀 分 布 随机 数组 ， 并 寻访 1) 数 
组 的 第 三 个 元 素 ; 2) 数组 的 1，2，5 个 元 素 组 成 的 子 数 组 ; 
3) 前 三 个 元 素 组 成 的 子 数组 ; 4) 寻 访 除 前 三 个 元 素 外 的 其 它 
元 素 ; 5) 前 三 个 元 素 倒 、 6) 由 大 于 0.5 的 元 
素 构 成 的 子 数组 ; 7) 将 元 素 的 第 {1，4 个 元 素 赋 值 为 0。 


rand(state ,0) 
Xx=frand(1,5) 
a1=X(3) 
a2=Xx([1 2 5]) 
a3=X(1:3) 
a4=X(3:end) 
a5=X(3:-1:1) 
a6=Xx(find(Xx>0.5)) 
X([1 4j)=[0 9] 





A=l123;456;789|; 


A(5,5)=11 1 A-[ 

A(:,6)=<<< 0 0 0 0 0 222 
AA=A(:,[1:6,1:6]) “0 
AA1=repmat(A,1,2) IEE 示 ， ， ， ， ， 
AA2=repmat(A,2,2) 0 0 0 0 11 222 


B=ones(2,6);AB_R=[A;B] ] 
AB _C=[A,B(:,1:5)7] 
s=[13689111416]; 

Al(S)=0O:A 


这 辑 1 标识 





邮 辑 过 辑 、 天 系 运 企 














any 一 吉 问 量 或 数组 的 列 癌 量 任意 元 际 不 为 0 则 返回 真 
all 一 奋 问 量 或 数组 的 列 癌 量 所 有 元 素 不 为 0 则 返回 其 
find 一 寻找 非 零 元 素 坐 标 

ULEUeIwtke) a=reshape([1:9],3,3); 

al(:,3)=Zeros(5,1) a=1./a 

a1=all(a(:,1)<10) f1=find(al) 

a2=all(a>3) 


f2=find(abs(a)>0.20&abs(a)<0.40) 
a3=any(al(:,1)>10) 


a4=any(a>10) 





2 了 22 
4=| 1 03 8 94 
| 了 3 7 2 


A=[265D5222) 1 63894) 13 457 21]，; 
alL=A(7) ，; 
下 这 天 ] 这 臣 二 直 7 已 工 ) 
a2=max (max(A) ) ; 
玉生 荆门 生 人 
7 引 2 ) 
[R V]=fEinaq(A==a2) ， 
fprIntE( 
rRyrV) 
L=aps (A) >10; 
X=A ( 工 ) 


和 天 阵 的 特殊 操作 





人 1) 矩阵 变 维 
。 freshape(X,M,N) 
2) 矩 阵 变 向 
。 rot90 (A) ，rot90 (A,K ) 
eelleigy ieieli 人 所 eli 
3) 短 阵 的 抽取 
。 diag(X, 八 ) 
。 tril(X)，tril( 久 ,人 )，triu(X)，triu( 义 , 八 ) 





抵 阵 的 基本 ,性质 





Size - Size of array. 

length Te 区 es ete] 

el Ne ee Eesti[eliTS 

numel - Number of elements. 

Isempty - True for empty array. 

>elll - True if arrays are numerically equal. 


Isequalwithequalnans -True ifarrays are numerically equal. 





Matlab 短 阵 巴 数 





天 阵 分 析 
norm - 天 阵 或 向 量 范 数 ; 
rank -和 珑 阵 秩 ; 
e[:] -上 短 阵 的 行列 式 
线性 方程 组 求解 
\ 和 / - 线性 方程 来 解 
inv -天 阵 来 间 
Con -天 阵 条 件数 
过 LU 分 解 
特征 值 和 特征 向 时 
eig -矩阵 的 特征 值 和 特征 向 量 
Ye - SVD 分 解 





左 除 和 右 除 的 比较 








用 以 下 程序 生成 条 件数 很 大 的 方程 

zanadn (' state' ,0) ; 

A=galLLezry ('zandsvd' ,100,2e13,2) ; s 产 生 条 件数 为 2e13 的 100 阶 窍 阵 
X=ones (100 ,1I) ; gs 生成 真 解 

DP=RAx*x; s 用 RAR 和 x 生 成 b 

比较 左 除法 和 求 逆 法 的 求解 所 用 的 时 间 和 相对 残 差 


tiC % 求 逆 法 

Xi=inv(A)”*b; 

ee 用 =ETeldiTI 人 人 人 丙 [ETeliiiifgabae)7dediile) 

Te % 左 除法 

《ee 

le 呈 lee 有 -ie 局 Wedle 人 eejLs ELeldilgesbhKeaej1aieliiiite) 
结 采 : 


ti = 0.1950, eri = 0.0539, rel =0.0037 
td =0.0156, erd =0.0779, rel =2.8535e-015 


本 章 开 始 问 题 的 求解 


求解 程 厅 : 





A=[0.07 0.18 0.15 0.24; 
0.04 0.24 0.10 0.65; 
0.54 0.42 0.54 0.10; 
0.35 0.16 0.21 0.01]; 
B =[0.15”*70; 0.25 ”70; 0.40 70; 0.20”70]; 
X=A\B 





反应 天 阵 的 秩 与 独立 反应 数 








独立 反应 : 在 一 个 存在 多 个 反应 的 复杂 反应 体系 中 ， 某 些 反 应 的 化 学 计量 方程 
可 以 由 其 它 反应 的 化 学 计量 方程 线性 组 合 得 到 。 如 果 m 个 同时 发 生 的 反应 中 ， 若 每 
一 个 反应 的 计量 方程 都 不 能 有 其 它 反 应 的 计量 方程 的 线性 组 合 得 到 ， 则 称 这 m 个 反 
应 是 相互 独立 的 。 


反应 扼 阵 : 


NeNNeN ReNRNLeX 





b 王 rank(A) 
解 得 ，b 二 3， 因 此 只 
办 有 三 个 反应 是 独立 的 








线性 方程 组 的 迭代 解法 


雅 可 比 友 代 : 


ie 





。1< w <2 用 于 加 速 某 收敛 的 迭代 过 程 ， 而 取 0< wo <1 用 于 非 收 


敛 欠 代 过 程 使 其 收敛 。 
“最 佳 松弛 因子 的 选取 ， 一 般 需 在 计算 过 程 中 搜索 寻 优 





线性 方程 组 的 迭代 解法 


在 偏 微分 方程 差分 解法 中 出 现 如 下 典型 方程 组 ， 试 
用 分 别 用 雅 可 比 和 迭 代 ， 斯 - 赛 德尔 迁 代 和 SOR 法 求 


ft EREIIXEIESTT 











一 
function Cha3demo7 
A=[-4 10000;)1L1 -41LI000;) 01 -41 0 0;0 0 工 
-410;0001 -41;0000 1 -4]:; 
D=diag (qiag (A) ) ; U=-t 上 ziu(A,，IL) ; 工 =- 七 ZI (A,， 一 工 ) ; 
XO=O*ones (6) ' :; 
pb=[=27 一 5 一 5 一 5 一 15 一 上 5 ， 
当 5JacobLan 工 ezatILon 
BJ=DAN (L+U) ; gJ=DN\Pb; Y=BJxx0O+gJ; 
mn= 荆 ; 
whiIle norm (Y-XO0) >=e 一 6 
X0=yY; 
Y=BJx*xO+Gg9J，; 
mn=Dm+ 二 + : 
end 
QisP('TIThe soLution of Jacobian It 上 eratILIon are ') 
fEPzintEf('N\nn=s3dN\n' yn) 
disP('y='),disp(y') 


ssGauSss 一 SelLadel 
BG= (D-L)N\U;gG= (D-L) \Pb;yYy=BGx*x0+GgG:; 
mn= 荆 ; 
noLrm (Y 一 X0) >=1Le 一 8 
X0=Y; 
Y=BG*xO+GgG ; 
mn=Dm 二 + : 


QIsP (人 

) 
fEPintE( 7 卫 ) 
QisP ) ,disP(Yy ) 


当 SOR 
omiga=1.08;BS= (D-omigaxL)N\((1IL 一 
omiga) xD+omigaxU) ;gGS=omigax ((D--omigaxL)N\P) ; 


Y=BSx*Xx0O+GG ; 
mn= 工 ; 
norm (Y 一 X0) >=1Le 一 8 

Xx0=Y:; 

Y=BSx*Xx0O+G9S ; 

mn=mn 二 LI :; 
disP( ) 
EPZIntE ( 7 卫 ) 


disP ) vdisP(y ) 


第 四 讲 
播 值 、 拟 合 与 数值 微分 和 积分 


化 工学 院 软 件 应 用 教科 组 
2006-10 


本 章 知 识 要 点 





。 插值 方法 (interp,spline) 

。 拟 合 方法 (polyfit,csaps) 

。 数值 微分 人 

。 数值 积分 (quad, quadl, fnint) 


elelNe[= ie 


插值 、 拟 合 、 数 值 微 分 、 数 值 积分 
在 化 工 计 算 中 的 作用 

。 表 格式 物性 数据 的 内 插 

。 离 散 实 验 数据 点 的 处 理 

。 状 态 方 程 计 算 流 体 的 烩 和 炳 

"微分 法 反应 动力 学 方程 拟 合 

。 等 温 活 塞 流 反应 器 的 设计 计算 
“微观 离 析 反 应 器 的 计算 





插值 向 介 


插值 的 ; 问题 可 以 描述 为 : 已 知 n+1 个 数 对 {xi, f(xi)}， 其 
中 ij=0,1...n， (Xi 互 不 相同 ， 称 之 为 节点 ) ， 求 取 函 数 
Jo 

当 {xi, f(xi)} 有 相当 的 精确 度 ， 但 它们 的 阴 数 关系 难以 确定 或 
难以 计算 时 ， 则 可 利用 这 些 数据 点 来 构造 一 个 较 简 单 的 函 
数 来 近似 表达 原 函 数 关系 。 

根据 通 近 函数 的 不 同 ， 第 见 的 插值 方法 : 
“Lagrange 多 项 式 插 值 〈 线 性 插值 ) 

“分 段 插值 


.有 理 式 插 信 





density~temperature 














180 


1 1 1 
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viscosity~temperature 


1 
470 
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Matlab 的 插值 (Interpolation ) 函数 





一 维 插值 使 用 FFT 方 法 的 一 “| interpft 
维 插值 
快速 一 维 插值 分 自 = 次 Hommite 本 


二 维 插值 一 


三 维 插值 er 





一 维 播 值 interp1 


调用 格 却 : 
下 = interp1(Xx,y,xi) 已 知 数据 向 量 


X,y ) ， 计 鼻 并 返回 在 插值 向 量 


Xi 处 的 函数 值 
yi=interp1(x,y,Xxi，method  ) 
yi=interp1(x,y,Xxi，method ， 
extrap ) 

Imethod' 用 于 指定 播 值 算 法 ， 其 
值 可 以 走 : 
nearest 一 一 节 近 插值 
inear 一 线性 插值 (默认 值 ) 








spline 一 一 分 段 三 次 样 条 插值 
pchip' 分 段 三 次 Hermite 插 值 


ee [ee 二 eielijje 几 局 





注意 : 
vv 向量 x 为 单调 。 若 y 为 矩阵 ， 
则 对 y 的 每 一 列 进行 插值 
v 向 量 xi 中 有 元 素 不 在 X 的 范 
转 内 ， 则 对 应 yi 值 为 NaN 
Y extrap 用 于 指定 当 向 量 Xi 
中 有 元 素 不 在 X 的 范围 内 
时 ， 采 用 "method 所 指定 的 
插值 算法 进行 外 播 计 算 与 之 
对 应 的 yi 值 





一 维 播 值 方法 比较 
已 知 X=[0:10]，y=[0 ”0.8415 “0.9093 ”0.1411 -0.7568 - 
0.9589 -0.2794 0.6570 0.9894 0.4121 -0.5440];(y = 
Sinx)， 比 较 一 维 线性 、 线 性 最 近 、 立 方 和 三 次 样 条 播 值 所 得 Xi = 
0.0.15,0.30.,0.45,.,10 处 的 值 yi。 如 果 初 始 数据 点 为 X= 
0,2,4,...,10，y = SinX， 以 上 方法 播 值 效 果 。 
Xx=0:1:10; 
y =[0 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.657 0.9894 0.4121 -0.5440]; 


ellt 全 全 we 有 Leegeli 
el et(CEILURLE 


Xi=0:0.15:10; 

yi=interp1(X,y,Xi); 
plot(xi,yi,r+),text(0.7028,0.4649,' 线 性 插值 rightarrow'") 
yi2=interp1(Xx,y,Xxi,nearst); 
plot(xi,yi2,"c*),text(3.537,0.1374,Neftarrow 了 最 近 插 值 ') 
yi3=interp1(X,y,Xi,cubic'); 
plot(xiyi3,,md'),text(2.408,0.8333,Neftarrow 三 次 插值 ') 
yi4=interp1(X,y,Xxi,spline ); 

plot(xi,yi4,"kh'),text(4.62,0.8158,' 三 次 样 条 插值 Yightarrow') 


初始 数据 对 于 播 值 的 影响 








Xx=0:2:10;y=Sin(X); 

el 全 ele 双 及 ieliegell 

ydelelt(CL LE 

Xi=0:0.15:10; 

yi=interp1(X,y,Xi); 
plot(xi,yir+),text(0.5876,0.2537,Neftarrow 线 性 插值 ') 
yi2=interp1(x,y,Xxi,nearst); 
plot(xi,yi2,"c*"),text(6.947,-0.258,Nleftarrow 也 近 插值 ') 
yiI3=interp1(X,y,Xxi, pchip ); 

plot(xiyi3mdy text(2.408,0.8333 ,Aleftarrow 三 次 插值 
yi4=interp1(Xx,y,Xi,Sspline ); 
plot(xi,yi4,"kh),text(1.601,1.138,Neftarrow 三 次 样 条 播 
值 ) 








Se eeile 


Spline() 的 调用 格式 为 : 

yi=spline(x,y,xi) “此 函数 等 同 于 yi=interp1(X,y,Xi， 
Spline ) 

pp=Sspline(x,y) 返回 三 次 样 条 插值 的 分 段 多 项 式 
形式 的 向 量 

spline 函 数 可 以 保证 插值 函 数 的 三 阶 手 数 连 续 


pchip() 的 册 用 格式 为 : 

yi=pchip(x,y,xi) “此 函数 等 同 于 yi=interp1(X,y,Xi， 
pchip ) 

pp=pchip(x,y) ”返回 三 次 样 条 插值 的 分 段 多 项 式 形 
式 的 向 重 





pchip 与 spline 的 区 别 


pchip 与 spline 的 区 别 : 三 次 样 条 在 相 邻 的 节点 上 并 不 保证 单 
调 性 ; 而 Hermite 分 段 三 次 样 条 则 可 保证 播 值 的 局 部 单调 性 


利用 点 (X=Sin(kK 工 /6),y=cos(k 关 /6)， 其 中 k=[0 1 2 3] 来 逼近 单 
位 圆 的 前 四 分 之 一 圆周 。 比 较 pchip 与 spline 的 差别 。 


t=linspace(0,pl2,4) 

Xx=Ccos(t);y=Sin(t); 

XX=|inspace(0,1,40); 

plot(X,y,S ,xx,[pchip(X,y,Xxj;Sspline(X,y,XX)]) 
eegeli 本 :MES=ielgl 

[elelle@lieliiI etcitcieeailieselllilis 








csape() 的 调用 格式 为 : 
pp=csape(x,y) 
pp=cSsape(X,y,Conds) 


conds 可 为 以 下 字符 串 : 

complete' 指定 端点 处 一 阶 导 数 

eeliiel 指定 端点 ER。 

'periodic 右 端 点 处 一 、 二 AS 
not-a-knot' 二 个 和 倒数 第 节点 处 三 阶 导数 连续 
variational 自 然 边 界 条 件 ， 即 六 点 全 从 王 数 为 二 


除 csape 外 ，Matlab 的 样 条 工具 箱 还 提供 了 其 它 样 条 插值 函 
数 ， 如 csapi，spapi 等 


CSsape 的 使 用 





设 f(X) 为 区 间 [0,3] 上 的 函数 ， 剖 分 节点 为 Xi=[0 1 2 3]; 节点 上 
的 函数 值 为 yi=[0 0.5 2 1.5]， 左 右 端点 处 的 一 阶 导 数值 分 别 
为 0.2 和 -1。 试 来 区 间 [0,3] 上 满足 上 述 条 件 的 三 次 样 条 插值 函 


程序 : 

Xx=[0 1 2 3]; 

y=|[0.2 0 0.5 2.0 1.5 -1]; 
6 二 ee 二 和 和/ 人 welillelisils 





几 个 有 用 的 函数 3 


1.[breaks coff k m nj]=uUnmkpp(pp) 命 令 可 以 看 到 样 
条 函数 的 具体 信息 
其 中 coff 是 一 个 矩阵 ， 其 第 i 行 是 第 i 个 三 次 多 项 式 
的 系数 ， 多 项 式 形式 如 下 : 





2.fnval 或 ppval 可 以 计算 样 条 函数 在 指定 点 的 函数 值 
3.fnder 和 fnint 可 分 别 计 算 样 条 函数 的 导数 和 积分 
4.fnplt 可 用 于 绘制 样 条 函数 的 图 形 


二 维 插值 : interp2 





调用 格 陈 : 


zl=interp2(Xx,y,Z,Xxiyl method ) 


Imethod' 算 法 属性 值 可 以 是 ; 





nearest 最 近 播 值 
由 inear 一 线性 插值 (默认 ) 
spline 一 一 三 次 样 条 播 值 (spline) 


cubic 一 立方 插值 





二 维 插值 函 数 的 使 用 WA 
号 一 组 分 度 系 数 的 “海底 深度 测量 数据 "， 由 以 下 一 段 程 邦 


randnfstate ,2); 

Xx=-5:5;y=-5:5;[X,Yj=meshgrid(Xx,y); 
Z=-500+1.2*exp(-((X-1).^A2+(Y-2).^2))-0.7*exp(-(exp(X+2).^2+(Y+1).^2)); 
Surf(X,Y,Z)，Vview(-25,25) 


试 由 插值 方式 绘制 海 反 形状 图 。 


Xi=linspace(-5,5,50); 
yi=linspace(-5,5,50); 
[YI=meshgrid(xi,yD); 
ZI=interp2(X,Y,Z,XLYL*cubic)) 
surf(XLYLZI) 

View(-25,25) 








拟 合 简介 


。 拟 合 和 插值 的 区 别 在 于 : 

1. 拟 合 时 ， 所 得 函数 不 需要 过 所 有 数值 点 

2. 插值 函 数 不 宜 外 推 ， 拟 合 函 数 在 某 些 情况 则 可 以 
。 拟 合 方法 中 最 第 用 的 是 最 小 二 乘 曲线 拟 合 


。 最 小 二 乘法 的 基本 思路 是 使 拟 合 因 变 量 y 在 给 定点 Xi 上 
使 残 差 平方 和 最 小 


用 于 拟 合 的 函数 可 以 是 多 样 的 ， 本 讲 只 介绍 多 项 式 函 
数 拟 合 和 样 条 函数 拟 合 





最 小 二 乘 多 项 式 拟 合 : polyfit 


调用 格式 如 下 : 

pP=polyfIt(X,y,m) 

[p,sj=polyfit(X,y,m) 
说 明 : 
国 输 入 参数 : (Xx,y) 为 已 知 数据 向 量 ，n 为 多 项 式 阶 数 ; 输出 
参数 p 为 拟 合 生成 的 多 项 式 的 系数 向 量 (长 度 为 n+T) ，S 为 
结构 参数 ， 供 函数 polyval() 调 用 以 获得 误差 估计 值 . 
国 邓 数 polyval() 第 第 与 polyfitO) 联 合 使 用 ， 其 调用 格式 为 : 
eelNA《2L(e 太 4 
国 [y,deltay]=polyval(p,xis), 返回 Xi 处 的 拟 合 函数 值 ，deltay 
是 50% 置 信 区 间 的 y 的 误差 。 


多 项 式 次 数 对 拟 合 效果 的 影响 








已 知 下 表 数 据 ， 用 polyfit 进 行 多 项 式 拟 合 





Xx=0.5:0.5:3;y=[1.75 2.45 3.81 4.80 8.00 8.60]; 
Xi=0.5:0.05:9; 
[a2,S2]=polyfit(x,y,2);[y2,delta2]=polyval(a2,xi,S2); 





plot(x,y,ro'Xxiy2,"m-),text(1.04,4.67, 二 次 多 项 式 拟 合 \eftarrow 妆 ) 

[elleEe]4 

[a4,S4]=polyfit(x,y,4);[y4,dqelta4]=polyval(a4,xi,S4); 

plot(x,y,ro'Xxiy4,"b-),hold on,text(1.563,6.775,' 四 次 多 项 式 拟 合 \leftarrow 妆 ) 
[a7,S7Z]=polyfit(x,y,7);[yZ,delta7]=polyval(a7,xiS7); 

plot(x,y, ro',xijy7,k-),hold on,text(1.931,9.532, 七 次 多 项 式 拟 合 \eftarrow'),hold off 


最 小 二 乘法 拟 合生 成 样 条 曲线 








函数 名 拟 合 准则 





。 拟 合 得 到 曲线 函数 sp 以 后 ， 可 利用 fnval(0) 计 算 任 
意 目 变量 下 的 函数 值 。 





函数 csaps() 的 用 法 CC 


局 


功能 : 平滑 生成 三 次 样 条 函数 ， 即 对 于 数据 (Xi,yi)， 所 求 的 三 次 样 条 
马 数 y=f(X) 满 足 : 





调用 格式 : sp=csaps(Xx,y,p) 
YyYS=CSapS(X,y,P,XX,W) 
输入 参数 : X,y “要 处 理 的 离散 数据 (Xi,yi) 
p 平滑 参数 ， 取 值 区 间 为 [0,1]。 当 p=0 时 ， 相 当 于 最 小 二 乘 直 线 拟 
合 ; 当 p=1 时 ， 相当 于 “自然 的 ”三 次 样 条 拟 合 
XX “用 于 指定 在 给 定点 XX 上 计算 其 三 次 样 条 函数 值 (yS) 
W ， 权 值 ( 权 重 ) ， 默 认为 1 
输出 参数 : Sp ， 拟 合 得 到 的 样 条 函数 
YS 在 给 定点 XX 上 的 三 马 数 值 





马 数 Spap2 () 的 用 法 





调用 格式 : sp=spap2(knots,k,X,y) 
sp=Spap2(knots,K,X,y,W) 


输入 参数 : knots 节点 序数 (knot sequence) 
k 样 条 函数 的 阶 次， 一 般 取 k=3， 有 时 取 k=4 
xy 要 处 理 的 离散 数据 (xi yi) 
W 权 值 (权重 ) ， 默 认为 1 


输出 参数 : sp 拟 合 得 到 的 样 条 函数 





马 数 Spaps () 的 用 法 


功能 : 平滑 生成 B 样 条 函数 





调用 格式 : sp=spaps(x,y,tol) 
[sp,ysj=spaps(xy,tolm,w) 


输入 参数 : X,y 要 处 理 的 离散 数据 (Xiyi) 
tol 光滑 时 的 允许 精度 
m 默认 值 是 2， 即 平滑 生成 三 次 B 样 条 曲线 
由 权 值 (权重 ) ， 默 认为 1 

输出 参数 : sp 拟 合 得 到 的 样 条 有 函数 


yS 在 X 上 经 平滑 处 理 的 B 样 条 函数 值 


马 数 csaps( 的 用 法 








采用 样 条 拟 合 函 数 csaps 进 行 拟 合 ， 比 较 p 的 取 值 对 于 拟 合 
效 采 的 影响 





X=0.5:0.5:3; 
y=[1.75 2.45 3.81 4.80 8.00 8.60] 
Xi=0.5:0.05:3; 
lot(X,y, ro 

4 对 于 微小 的 数据 扰动 ， 多 项 
c=[000:100:001;011;0.50.50.5]， 式 拟 合 通 第 比 样 条 函数 更 为 
| 敏感 
for i=0:0.25:1.0 

Oj=Cd)， 


yS(,:)=CSapS(X,y,XD); 
plot(Xxi,ys(j,2), color ,Cj),pause 
j=j+1 

ie 


数值 微分 





。 对 于 列表 型 函数 往往 需要 用 数值 方法 计算 函数 的 微分 
。 数值 微分 的 基本 方法 
1. 兰 分 
2. 利用 插值 ( 拟 合 ) 多 项 式 来 微分 
3. 利用 三 次 样 条 插值 ( 拟 合 ) 函数 求 微分 
。 数值 微分 可 以 放大 误差， 应 谨慎 使 用 





Matlab 数 值 微分 实现 方法 Ce 





有 限 差 分 法 : 用 差分 函数 diff() 近 似 计 算 导 数 

多 项 式 拟 合 方法 : 
人 

三 次 样 条 插值 方法 
离散 数据 ED 实 杀 拓 值 马 次 cs 加 罗 了 二 次 pp EDEZO 在 xi 的 导 效 人 
样 条 拟 合 方法 


误 族 数 失 上 ET 了 国 if 条 合 己 交 sp 国有 :2 交 pp ED 在 xi 的 导 效 侍 





画 数 diff 


对 于 向 量 X，diff(X) 表 示 了 [X(2)-X(1) X(3)-X(2) .… X(n)-X(n-1)] 
对 于 矩阵 X，diff(X) 表 示 了 [X(2:n,:) - X({:n-1.]] 
diff(x,n,dim) 得 到 短 阵 x 在 dim 维 上 的 n 阶 差 值 


>> diff((1:10).^2) 一 一 ww [135791113151719] 


x=[138;246] 
diff(x,1,1) diff(Xx,1,2) diff(X,2,2) diff(Xx,3,2) 
星 星 星 
1 1 -2 3 
- [2 5ez| Empty matrix: 2-by-0 
el [et (Ce 和 【9 
Jel egell 
利用 diff 画 数 求 sin(X) ww 
让 el ef 


eNet(STiLe4e4)740TE 
plot([1:h:10-hj,diffy [hh:10-hhj,diffyy, KK-o)) 


邓 数 gradient 


FX = gradient(F) 
Eee 
[Fx,Fy,Fz] = gradient(F ) 
[| = gradient(F ,) 

[..] = gradient(F,h1,h2…) 


下 

1，Fx=dF/dx，FY=dF/dy，Fz=dF/dz; 当 F 为 向 量 时 ，dF=gradient(F) 是 一 维 
梯度 

2 N2 二 是 步 长 ， 缺 省 值 为 1 


3， 如果 h1，h2，..….. 为 向 量 ,其 长 度 必 须 匹 配 F 的 维 数 





例题 





某 一 液 相 反应 浓度 随时 间 变 化 的 实验 数据 如 下 表 : 


TooTOT2T50TOD 





ce 15191377 1230|157 | 08 | 025 |0094 


试 t=0.1,0.4min 时 的 反应 速度 


x=[0 0.2 0.6 12 5 10]; 

C=[5.19 3.77 2.30 1.57 0.8 0.25 0.094]; 

ell 人 Geieo 有 ielegeil 

esizleStaeNUelttee) 

eg 二 ielzlifeje) 

ee 的 器 ee:L(eLGRR 

LTG 

fprintf(C The reaction rate at 0.1min and 0.4min are2%5.4f,%e5.4f... 
(LIDILS ea NAUEeL 的 区 eg 


数值 积分 





。 数值 积分 在 数值 计算 中 有 着 重要 作用 ,许多 数值 计算 问题 
可 以 转化 为 数值 积分 问题 ,如 常 微分 方程 初 值 问题 等 
。 通 常 可 用 逼近 多 项 式 Pn(X) 来 代替 被 积 函 数 f(X)， 计 算 积 


N 
分 






。 构造 数值 积分 的 方法 很 多 ， 主 要 有 Newton-Cotes 系 列 
数值 积分 法 、Gauss 积 分 法 和 Romberg 积 分 法 等 


数值 积分 





TOTTTE 
自 适应 Simpson 求 积 公式 〈 低 阶 ) 


用 


单 形 求 积 公 式 ， 速 度 决 ， 精 度 关 


弟 形 法 求 一 个 区 间 上 的 积分 曲线 


人 个 区 间 上 的 积分 曲线 ， 精 
度 仿 专 

fnint 利用 样 条 函数 求 不 定 积 分 ， 与 spline， 
ppval 配 合 使 用 ， 主 要 应 用 于 表格 “ 郴 数 ” 
积分 

















调用 格式 : ，z=trapz(y) 


> 用 梯形 来 积 方法 计算 y 的 积分 近似 值 。 

> 对 于 向 量 Yy，trapz(y) 返 回 y 的 积分 ; 

> 对 于 短 阵 Yy，trapz(y) 返 回 一 行 向 量 ， 向 量 中 的 各 元 素 为 珑 
阵 y 的 对 应 列 向 量 的 积分 值 ; 





目 运 应 Simpson 法 数值 积分 : quad 





调用 格式 : ”qdq=quad(@fun,a,b) 
e 忆 eLLELeLKCLLE2ReRe1 
qdq=quad(Qfun,a,b,tol,trace,p1,P2,………… ) 


输入 参数 : “fun ”被 积 了 数 。 在 定义 fun 时 ， 被 积 函 数 表 达 
式 必 须 是 向 量 形式 ， 即 表达 式 人 必须 使 用 点 
运算 符 (.*、./ 和 .A) 以 支持 向 量 
ab 。” 即 积分 限 [ab] 
el] 绝对 误差 限 ， 默 认 值 为 1.e-6 
p1,p2..…… 直接 传递 给 函数 fun 的 已 知 参 数 


输出 参数 : “qq 积分 结 采 


自 适 应 Lobatto 法 数值 积分 : quadl0 调用 格式 同 quad 


不 同 积分 函数 的 比较 








TIeiielEEGEEIelsiiliers 


、 eteliie 
术 和 分 [人 


nt=length (tb ; 


y=fun(t); 
比较 CU miISUm)， sc=cumsum(y)*d; 
Scf=sc 人 nt) 
trapz ? e[I[UeE 了 本 
ee 二 si| un,0,3*pi,[， 
eictel 的 积分 精 qd2=quadl(@fun,0.3*pi,[] 
、 EV=0.9008407878; 
大 ， 访 积分 个 的 ] 精 磺 err(1)=abs(scf-EV): 
刀 、 err(2)=abs(z-EV 六 10; 
解 为 err(3)=abs(qd-EV)*10; 
err(4)=abs(qd2-EV 六 10; 
0 9008407878 bar(ern),title( 不 同 积分 方法 比较 


weleldiuiiietulsg 

text(0.7,7e-4,' 珑 形 法 ',fontsize',20) 
text(1.5,5e-4,' 梯 形 法 ,fontsize',20) 
text(2.4,3e-4,,'Simpson 法 ,fontsize',20) 
text(3.4,1e-4,'Lobatto 法 ',fontsize',20) 


ITea LT 人 
y=exp(-0.5"t).*Ssin(t+pl6) 


广义 积分 


1. 奇 点 积分 2. 无 穷 积 分 


fun=inline(1./(sgqrt(X).*(exp(X)+1)); 可 先 选 取 一 个 有 限 的 积分 区 间 ， 如 
etEtelitIRD 


[0,100] 计 算 ; 在 选择 一 个 较 大 的 积 

quadl(fun,eps,1) 分 区 间 ， 如 [0 200] 计 算 ， 如 两 次 计 
算 结果 的 差 满足 一 定 的 精度 要 求 ， 
则 可 认为 此 值 即 为 无 穷 积 分 的 值 











多 重 数值 积分 Ce 





二 重 积 分 函数 : 


eolelitteitTR4LRBE 人 AL Are IT 

三 重 积 分 函数 : 
iteUBLRBERAAITIRIE 全 41DR4E 人 Ce] 几 IITIeieh 
说 明 : 

Matlab 只 能 处 理 积分 限 为 常数 的 多 重 积 分 ， 对 内 积分 上 下 限 为 
外 积分 变量 的 积分 问题 , 需 自行 编写 函数 求解 。 





样 条 函数 在 数值 积分 与 微分 中 的 应 用 


已 知 X=(0:0.1:1)*2*pijy=sin(X) ， 求 其 积分 和 导数 


X=(0:0.1:1 六 2*pijy=Sin(X); pp=spline(Xx,y); 
Int_pp=fnint(pp); 
e[=] 且 ee 呈 册 elifeje) 


9% 在 基础 区 间 上 ， 计 算 三 个 样 条 函数 与 理论 值 的 最 大 误差 

X=(0:0.01:1 六 2*pi; 
EM 人 leS(ejeifee 太 44430447) 
EMtzleSfeeLtil 本 ee 太 44 苹 人 weib44)) 
3 本 o3 呈 UEMNMEleSilejel(els ee 及 444)) 
9% 计算 y(X) 在 区 间 [1,2] 上 的 定 积分 


了 siilitsililteieltsieNAsjellii[s eeLL ee 全国 的 有 
Definitelntegral.byTheory=(1-cos(2))-(1-cos(1T)); 

9% 计算 dy(3)/dx 

BT NeNeseliils ae ee 

Derivative.byTheory=cos(3); 

Derivative.byDiference=(sin(3.01)-sin(3)X0.01; % 前 向 差分 近似 
[Biilaliteiiiltsreigz 四 BIN 

eltteje 同 Nielegell 

ent 本 ee 区 册 国有 ijellt(els ee 有 二 ielegeil 

legend(Cy(xX) ，S(x dy/dx)) 






反应 器 停留 时 间 分 布 的 混合 特性 多 


在 t- 0 的 时 刻 ， Ni 
入 一 定量 的 示 足 剂 ， 同 时 在 容器 出 口 处 测量 流出 物料 中 示 踪 剂 
浓度 随时 间 的 变化 ， 实验 数据 如 下 表 : 











EREIEIEIEIEIEIIEIEIICIEIET 
EREIEIENERIIIEGECSETTITIIG 


试 计算 流体 在 容器 中 的 平均 停留 时 间 以 及 扩散 准 数 。 


数学 模型 
平均 停留 时 间 多 扩散 特征 数 为 : 


LaieEeicEre site 

t = [0:20:200];C = [0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0];plotlt,C,o) 
sp1= splinelt,C); 

[elleEeil 。 
fnplt(sp1,b-); xlabel(Time (s)))，ylabel(C (kmoymA^3)iA10^31) 
LeleEe il 

% By spaps(): 

figure，plotlt,C,'o1) 

Sp2= Spaps(t,C,1); 

eegell 。 
fnplt(sp2,b-); xlabel(Time (S))，ylabel(C (kmolym^3)iA10^3) 
eegelii 

器 eGAAwiEILILeieESiellleilw Isieateilliltectrc 

% Integration 

to0 = Of = tend); 

IC = quadl(@Func:to,tf, 中 ,由 ,sp); 

ICct = quadl(@Func1,to, 蔷 ,由 [0,sp); 

ICt2 = quadl(@Func2,t0, 革 ,由 ,由 ,sp); 

tm1 = ICtUIC 

ss1 = |Ct2/IC - tm1^2 

DL2uL1 = ss1/tm1^2/2 


function y = Func(x,Sp)  %f= C 
和 三 由 X); 


LIwieliNAE 司 本 (x,SP) %f= Ct 

y = fnval(sp,X)…X; 

96 ------------------------------------------------------------------ 
function y = Func2(Xx,Sp) 9% f= Ct^2 

y = fnval(sp,X)…*X.^2; 


微分 法 进行 动力 学 数据 分 析 Ce 





反应 物 A 在 一 等 温 间歇 反应 器 中 发 生 反 广 为 瑟 2 人 

测量 得 到 反应 器 中 不 同时 间 下 反应 物 A 的 浓度 CA 如 

下 表 所 示 。 试 根据 表 中 数据 确定 其 反应 速率 方程 。 
EDIEIEIEIEEIIIEI 
al 





TIeielEGTLEielslilieisl 

9%6 动力 学 数据 

t=[020 40 60 120 180 300];CA=[10865321; 
% 用 最 小 二 乘 样 条 拟 合法 计算 微分 dCA/dt 

knots = 3:K= 3; 5% 三 次 B 符 条 

sp = Spap2(knots,K,tCA); 

sp = Spap2(newknt(sp), 人 ,ttCA); 

pp = fnder(sp); 9% 计算 B 样 条 函数 的 导 函 数 
dCAdt = fnval(pp 刷 ; % 计算 t 处 的 导 函 数值 

9%6 绘制 图 形 

ti = linspacelt(1)ttend),200); 

CAi = fnval(Ssp,th; 

plotlt,CA, ro tiCAib-)，xlabel(t)，ylabel(C_A) 
figure 

fnplttppP) 

Ge TIT(ejeR 

ettiIEeL7aelilR 

[21el=] [他 时 

EeeL@le 

26 线性 拟 合 

07el 

y=|og(CTA); 

Xx = |og(CA); 

p = polyft(Xy,1); 

ee 

n = p(1) 


第 第 五 寺 
第 微分 方程 数值 解 


化 工学 院 软件 应 用 教科 组 
2006-10 


本 章 知 识 要 点 





。 数 信 计算 
- 利 微 分 方程 初 值 问题 
- 币 微 分 方程 边 值 问题 


。 MATLAB 


- 徽 分 方程 求解 第 微分 方程 了 的 相关 困 数 
eele[72e) 
。bvp4c 





微分 方程 在 化 工 模型 中 的 应 用 





“间歇 反应 器 的 计算 

。 活 塞 流 反应 器 的 计算 

。 全 混流 反应 器 的 动态 模拟 

。 定 态 一 维 热 传导 问题 

逆流 壁 冷 式 国定 床 反应 器 一 维 模型 
固定 床 反应 器 的 分 散 模型 


Matlab 第 微分 方程 求解 问题 分 类 








初 值 问题 : 边 值 问题 : 

. 定 解 附加 条 件 在 自 变 量 > 在 自 变 量 两 端 均 给 定 附 加 
的 一 端 条 件 

. 一 般 形式 为 - 

。 初 值 问题 的 数值 解法 一 > 边 值 问题 可 能 有 解 、 也 可 
般 采 用 步 进 法 ， 如 能 无 解 ， 可 能 有 唯一 解 、 
Runge-Kutta 法 也 可 能 有 无 数 解 

>> 边 值 问题 有 3 种 基本 解法 
。 人 帮 加 法 
。 打靶 法 


。 松弛 法 


Matlab 求 解 稍微 分 方程 初 值 问题 方法 





将 待 求 解 转 化 为 标准 形式 ， 并 "翻译 ” 
成 Matlab 可 以 理解 的 语言 ， 即 编写 
odefile 文 件 


选择 合适 的 解 工 指令 求解 问题 


根据 求解 问题 的 要 求 ， 设 置 解 算 指令 
的 调用 格式 


二 国 
ode45 通 4-5 阶 法 解 ODE odeset 创建 、 更 改 ODE 选 项 的 
往 设置 


通 变 阶 法 解 ODE | 读 取 ODE 选 项 的 设置 
解 适度 刚性 ODE ODE 的 输出 时 间 序 列 图 
变 阶 法 解 刚性 ODE ODE 的 二 维 相 平面 图 
器 低 阶 法 解 刚性 ODE 机 和 x 间 图 
一 刚性 ODE 在 Mallab 信 人 窗 显示 结 








ele[ls 让 = 


CG 六 的 odefile 实 际 上 是 一 个 Matlab 函数 文件 ， 一 般 作 
整个 求解 程序 的 一 个 子 函 数 ， 表 示 0ode 求 解 问题 
C wab 了 odefile 的 模板 ， 采 用 type odefile 命 令 显 
示 其 详细 内 容 ， 然 后 将 其 复制 到 脚本 编辑 窗口 ， 在 合 
和 适 的 位 置 填 入 所 需 内 容 
G 一 般 而 言 ， 对 于 程 厅 通 用 性 要 求 不 高 的 场合 ， 只 需 将 
原 有 模型 写成 标准 形式 ， 然 后 “翻译 "成 Matlab 语 言 即 可 





odefile 的 编写 规定 


1.，ode 文 件 的 最 简单 格式 必须 有 一 个 自 变 量 t 和 
函数 y 作 为 输入 变量 ， 一 个 y 的 导 函 数 作为 输 
出 变量 。 其 中 自 变 量 t 不 论 在 ode 文 件 中 是 否 
使 用 都 必须 作为 第 一 输入 变量 ，y 则 必须 作为 
第 二 输入 变量 ， 位 置 不 能 颠倒 。 

2.， 可 以 向 ode 文 件 中 传递 参数 ， 数 目 不 受 限制 


odefile 的 编写 


求解 初 值 问题 : 





站 {f=fun(X,y) 


fy-2”X/y; 


初 值 问题 : [md (O<x<]) Enction f-fun(xy) 


{=y + Y^2; 














常 微分 方程 组 与 单个 常 微分 方程 求解 方法 相同 ， 
只 需 在 编写 odefile 时 将 整个 方程 组 作为 一 个 向 量 


function f=fun(X,y) 

dy1dx = 0.04*(1-y(1))-(1-y(2))."y(1)+0.0001*(1-y(2)).^2; 
dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; 
[eeheN2eh4 


高 阶 微分 方程 odefile 的 编写 
(7 0D)y =< cos2 卫 全 


本 例 的 难度 : 





方程 系数 非 线 性 方程 高 阶 ， 非 标准 形式 
可 在 odefile 中 定义 方程 变形 : 令 y1=y; y2=Vy' 
则 原 方 程 等 价 于 : 


function f=fun(t,y) 


a=-exp(-t)+Cos(2*pi*t)*exp(-2*H): 4 
b=cos(2*pi*t); | 区 
人 =[y(2) 


-a yc) <-b y(T)+exp(ttb D]; 





解 算 指 令 的 使 用 方法 
调用 格式 : 





1， [TY]j=ode45(@fun, TSPAN,Y0O) 

6 忆 eie[icI(CIiIRIES NEEORejeiielit 
ecI(CIULRNIES IN 了 人 elelieli 0 国 二 2 和 
本 EeelscI(CLTLIRISI IN 了 人 NeeilleliE 汪 国 二 729 交 ) 


说 明 : 

尝 输出 变量 T 为 返回 时 间 列 向 量 ; 解 矩 阵 Y 的 每 一 行 对 应 于 T 的 一 个 元 素 ， 列 数 与 求解 
变量 数 相等 。 

学 @Qfun 为 函数 多 柄 ， 为 根据 竺 求解 的 ODDE 方 程 所 编写 的 ode 文 件 (odefile ) ; 

必 TSPAN = [TO TFINAL] 是 微分 系统 Y = Flty) 的 积分 区 间 ; Y0 为 初始 条 件 

学 0ptions 用 于 设置 一 些 可 选 的 参数 值 ， 缺 省 时 ， 相 对 于 第 一 种 调用 格式 。options 中 
可 以 设置 的 参数 参见 0deset 


学 P1，P2，... 的 作用 是 传递 附加 参数 P1，P2，..: 到 ode 文 件 。 当 options 缺 省 时 ， 应 
在 相应 位 置 保留 []， 以 便 正 确 传 递 参 数 。 





常 微分 方程 初 值 问 题解 算 指令 比较 


odel15Ss 基于 数值 差分 的 可 变 阶 方法 (BDFs， 低 
Gear ) 


二 阶 改进 的 Rosenbrock 法 


站 
蜂 
使 用 樟 形 规则 四 


ode23tb TR-BDF2 〈 隐 式 Runge-Kutta 法 ) 





ode 解 算 指 令 的 选择 (1) 





1 .根据 常 微 分 方程 要 求 的 求解 精度 与 速度 要 求 


求解 初 值 问题 : RE (0<x<]) 


比较 ode45 和 ode23 的 求解 精度 和 速度 





ode45 和 ode23 的 比较 





LeieliEGtcieliiley 

人 Geliljel21SiellEeiasIsiI ESeleitLiLie 上 iielilEeiels2is=1ileEeiel=72x6 琶 :ielNs] 
elduitctalleliie 

y0=1 ; 

tic,[x1,y1]=ode45(fun,[0,1],y0)it ode45=toc 
tic,[x2,y2]=ode23(fun,[0,1],y0):tode23=toc 

plot(X1,y1,b- ,Xx2,y2,,m-.)， 

xlabel(Xx'),ylabel(y )， 

[lleielB NBS NLeieciiellNLelN SN 


disp(CComparative Results at Xx=1 2 ); 
fprintf( AnODE45N\ttt y=%.8fnODE23Nttt y=%.8fnPrecisive Result .…. 
=9%.8fn ,y1(end),y2(end),1.7320508) 


function f=fun(X,y) 
{f=y-2*X/y; 





ode 解 算 指令 的 选择 (2) 
2. 根 据 第 微分 方程 组 是 否 为 刚性 方程 


如 果 系 数 答 阵 A 的 特征 值 连 乘 积 小 于 零 ， 
[ww 且 绝 对 值 最 大 和 最 小 的 特征 值 之 比 (刚性 
比 ) 很 大 ， 则 称 此 类 方程 为 刚性 方程 


[0 Eee， 由 | ,性 比 : 100/0.01 = 10000 


。 刚 性 方程 在 化 学 工程 和 自动 控制 领域 的 模型 中 比较 利 
见 。 








ode 解 算 指 令 的 选择 (2) 


。 刚 | 性 方程 的 物理 意义 : 第 微分 方程 组 所 描述 的 物理 化 学 变 
化 过 程 中 包含 了 多 个 子 过 程 ， 其 变化 速度 相差 非常 大 的 数 
量 级 ， 于 是 常 微分 方程 组 含有 快 变 和 慢 变 分 量 。 

“党 微分 方程 组 数值 积分 的 稳定 步 长 受 模 值 最 大 的 特征 值 控 
制 ， 即 受 快 变量 分 量 约束 ， 特 征 值 大 则 允许 步 长 小 ; 而 过 
程 趋 于 稳定 的 时 间 又 由 模 值 最 小 的 特征 值 控 制 ， 特 征 值 小 
则 积分 到 稳定 的 时 间 则 长 。 

。Matlab 提 供 了 不 同 种 类 的 刚性 方程 求解 指令 : ode15s 
ode23s ode23t ode23tb， 可 根据 实际 情况 选用 





刚性 第 微分 方程 组 求解 








LeliEeitzteielsiiile 

figure 

ode23s(@fun,[0,100],[0;1]) 

elie 区 eli 

ode45(fun,[0,100],[0;]) 
人 
function f=fun(X,y) 

dy1dx = 0.04*(1-y(1))-(1-y(2))."y(1)+0.0001*(1-y(2)).^2; 
dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; 

[eehe2 oh 


解 算 指 令 的 输出 控制 


1. [TY]= ode45(@fun, TSPAN,Y0,options,P1,P2,..) 

2. sol=ode45(@fun,TSPAN,Y0) 将 解 输出 指 结构 体 sol 中 ; SXINT = 
deval(SOL,XINT) 用 于 计算 解 在 XINT 处 的 值 ，XINT 必 须 位 于 区 间 [SOL.x(1) 
司 @] 酌 《(=ilej 着 

3. 在 无 输出 变量 时 ， 将 调用 默认 的 odeplot 输 出 解 的 图 形 。 

4. 除了 以 odeplot 形 式 输出 外 ， 还 可 以 以 odephas2， 和 odephas3 的 形式 输出 解 
向 量 的 二 维和 三 维 相 平 面 图 。 

5. 采用 以 下 语 名 options=odeset(outputfcn',odephas2) 可 以 将 输出 方法 改变 为 
平面 输出 

6. odeprint 输 出 来 解 过 程 每 一 步 的 解 





解 工 指令 的 options 选 项 


~ GO OO 人 








.RelTol - 相对 误差 ， 它 应 用 于 解 向 量 的 所 有 分 量 。 在 每 一 步 积 分 过 程 中 ， 第 i 


个 分 量 误 差 e(i) 满 足 : el()<=max(RelTol*abs(y(i),AbsTol(i))。 


. AbsTol - 绝对 误差 ， 若 是 实数 ， 则 应 用 于 解 向 量 的 所 有 分 量 ， 若 是 向 量 ， 则 


它 的 每 一 个 元 素 应 用 于 对 应 位 置 解 向 量 元 素 。 


. OutputFcn - 可 调用 的 输出 函数 名 。 每 一 步 计 算 完 后 ， 这 个 函数 将 被 调用 给 


出 结果 ， 可 以 选择 的 值 为 : odeplot，odephas2，odephas3，odeprint。 


.OutputSel - 输出 厅 列 选择 。 指 定 解 向 量 的 哪个 分 量 被 传递 给 DutputFcn。 
. MaxSetp - 步 长 上 界 ， 缺 省 值 为 来 解 区 间 的 1/10。 
.InitialStep - 初始 步 长 ， 缺 省 时 自动 设置 。 

. 采用 odeset 改 变 原 有 选项 的 值 





二 2 \ 多 人 一 三 、 Ai 
间 区 反应 器 中 串联 -平行 复杂 反应 系统 


在 间 歌 反应 器 中 进行 液 相 反应 制备 产物 B， 反 应 网 络 如 图 5-1 所 示 。 反 
应 可 在 180~260C 的 温度 范围 内 进行 ， 反 应 物 X 大 量 过 剩 ， 而 C, D 和 
E 为 副 产 物 。 各 反应 均 为 一 级 动力 学 关系 : T= -kC， 式 中 睫 凤 
已 知 参 数 : ku; = 5.78052 x 1010，ko> = 3.92317 x 1012，koa3 = 
1.64254 x 104，ku4 = 6.264 x 108，Eai=124670，Ea2=150386， 
Eas=77954，Ea4=111528。 初 始 浓度 : CA=1kmolym3， 其 余 物 质 浓 度 为 0。 
已 知 是 产物 B 收 率 最 大 的 最 优 反应 温度 为 224.6 

试 计 算 1) 在 最 优 反 应 温度 下 各 组 分 浓度 随时 间 的 动态 变化 ; 2) 最 优 反 应 时 
间 ; 3) 和 输出 产物 D 对 反应 物 浓 度 A 的 关系 图 。 














LTwITeL Eee ie7 





T = 224.6 + 273.15; R = 8.31434; k0 = [5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8]; 


Ea = [124670 150386 77954 111528]; 
% Initial concentration:C0(i), kmol/m^3 
co=[10000]itspan=[01e4]; 


opt=odeset( reltol ,1e-4,outputfcn , odephas2 ,outputsel ,[1;4]) 


[tcC] = ode45(@MassEdquations, tspan, C0,opt,k0,Ea;,R,T) 
?%e plot 

plot(t'C(1) rtC(2)， ktC(:3) bb- tiC(:4)，k--); 
xlabel(Time (S) ); 

ylabel(Concentration (kmolMm^3) ); 
legend(A','B',C';D) 

CBmax = max(C(:,2)); % CBmax 

yBmax = CBmax/C0(1)  ?eyBmax: 

index = find(C(:,2)==CBmax); 

t_opt = tindexX) % t_ opt: the optimum batch time, s 
96 ------------------------------------------------------------------ 
function dCdt = MassEquations(t,C,k0,Ea,R,T) 

2% Reaction rate constants, 1/s 

k = k0.*exp(-Ea/(R*T)); k(5) = 2.16667E-04; 

2% Reaction rates, kmoles/m3 s 


rA = -(k(1)+k(2))*C(1); rB = k(1)*C(1)-k(3) C(2);rC = k(2)C(1)-k(4)*C(3);rD = k(3)*C(2)-k(5) C(4); 


rE = k(4)*C(3)+k(5)*C(4); 
2% Mass balances 
dcCdt = [rA;i rB;i rC;i rD;TrE]; 


Matlab 求 解 边 值 问题 方法 : bvp4c 函 数 





1. 把 待 解 的 问题 转化 为 标准 边 值 问 是 本 本 生 

2. 因为 边 值 问题 可 以 多 解 ， 所 以 需要 为 期 望 解 指定 一 个 初 
始 猜 测 解 。 该 猜测 解 网 (Mesh ) 包括 区 间 [a, b] 内 的 一 组 
网 点 (Mesh points ) 和 网 点 上 的 解 S(X) 

3. 根据 原 微分 方程 构造 残 差 巴 数 区 芝 JOB 

4. 利用 原 微 分 方程 和 边界 条 件 ， 人 和 借助 迭代 不 断 产 生 新 的 
S(X)， 使 残 差 不 断 减 小 ， 从 而 获得 满足 精度 要 求 的 解 


Matlab 求 解 边 值 问题 的 基本 指令 





selliiaexeilili 几 | Penmeeis ) 
生成 byp4c 调 用 指令 所 必须 的 " 解 猜测 网 


ie ee GoesiIiReiewilRESel LaoelieI ee 和 
给 出 微分 方程 边 值 问题 的 近似 解 


Sxint = deval (sol，Xxint ) 
计算 微分 方程 积分 区 间 内 任何 一 点 的 解 值 





初始 解 生成 函数 : bvpinit() 


solinit = bvpinit 〈《Xxv,parameters ) 


> X 指 定 边界 区 间 [a,b] 上 的 初始 网 络 ， 通 常 是 等 距 排 列 的 〈1 x M ) 一 维 数组 。 
注意 : 使 X(1)=a，x(end)=b; 格 点 要 单调 排列 。 
> V 是 对 解 的 初始 猜测 
> Solinit〈 可 以 取 别 的 任意 名 ) 是 " 解 猜测 网 (Mesh ) ”。 它 是 一 个 结构 体 ， 
带 如 下 两 个 域 : 
依 Solinit.x 是 表示 初始 网 格 有 序 节 点 的 (1x M) 一 维 数组 ， 并 且 
solinit.Xx(1) 一 定 有 是 a，Ssolinit.x(end) 一 定 是 pb。M 不 宜 取 得 太 大 ，10 数 
量 级 左右 即 可 。 
人 Solinit.y 是 表示 网 点 上 微分 方程 解 的 猜测 值 的 (Nx M ) 二 维 数组 。 
solinit.y(: 让 表示 节点 Solinit.Xx( 处 的 解 的 猜测 值 。 





eehe2iei 人 eleL>iilieieiliEE el eelieIi ee 羽生 的 


输入 参数 : 
> odefun 是 计算 导数 的 M 函 数 文 件 。 该 函数 的 基本 形式 为 : dydx = 


odefun(xy,parameters,p1,p2，...)， 在 此 ， 自 变量 x 是 标量 ，yY，dydx 
是 列 向 量 。 其 基本 形式 与 初 值 问 题 的 odefile 形 式 相 同 
bcfun 是 计算 边界 条 件 下 残 数 的 M 函 数 文 件 。 其 基本 形式 为 : res = 
bcfun(ya,yb,parameters,p1,p2,.………… )， 文 件 输入 宗 量 ya,yb 是 边界 条 件 
列 向 量 ,分 别 代 表 y 在 a 和 b 处 的 值 。res 是 边界 条 件 满 足 时 的 残 数 列 向 量 。 
注意 : 例如 odefun 函 数 的 输入 宗 量 中 包含 若干 未知" 和 "已 知 " 参 数 ， 那 
么 不 管 在 边界 条 件 计 算 中 是 否 用 到 ， 它 们 都 应 作为 bcfun 的 输入 宗 量 。 
合 入 宗 量 options 是 用 来 改变 bvp4c 算 法 的 控制 参数 的 。 在 最 基本 用 法 
中 ， 它 可 以 缺 省 ， 此 时 一 般 可 以 获得 比较 满意 的 边 值 问 题解 。 如 需 更 改 
可 采用 bvpset 函 数 ， 使 用 方法 同 odeset 函 数 。 
输入 宗 量 p1，p2 等 表示 和 硕 望 向 被 解 微分 方程 传递 的 已 知 参 数 。 如 果 无 
须 向 微分 方程 传递 参数 ， 它 们 可 以 缺 省 。 





边 值 问题 求解 指令 ; bvp4c () B 号 


输出 参数 : 

输出 变量 Sol 是 一 个 结构 体 

售 Sol.X 是 指令 bvp4c 所 采用 的 网 格 节点 ; 

售 Sol.y 是 y(X) 在 Sol.Xx 网 点 上 的 近似 解 值 ; 

售 Sol.yp 是 Y'(X) 在 Sol.X 网 点 上 的 近似 解 值 ; 

令 Sol.parameters 是 微分 方程 所 包含 的 未 知 参 数 的 近似 解 值 。 当 
被 解 微分 方程 包含 未 知 参 数 时 ， 该 域 存在 。 


边 值 问题 的 求解 到 








人 TeieliEGiiEtselsiliet 

。 寺 这 渴 区 关 作 | solinit=bvpinit(linspace(0,1,10),[1 0]); 
ie eeiew(CI@1B aiIRC11G31LRSielllilli 
X=[0:0.05:0.5]; 


原 方程 组 等 价 于 以 下 标准 形式 y=deval(sol,x); 
的 方程 组 : yP=[0 -0.41286057 -0.729740656.. 


-0.95385538 -1.08743325 -1.13181116]; 


XP=[0:0.1:0.5]; 

plot(XP,yP,，o ,X,y(1，, Dr-)， 

[elieleatiiMiesiESelILeNLLILITwsLESIellIiel 
Am 26 -一 


边界 条 件 为 dydx=[y(2)y(D)+10] 
26 





etelialeeaasi@iilTA《 枚 4e) 


Eeey》 bc=[ya(1)yb(1)] 


昌 
有 


边 值 问题 的 求解 








X = [0:0.12:1]; 
y = deval(sol,X); 
yP=[1 1.0743 1.1695 1.2869 1.4284... 
1.5965 1.7947 2.0274 2.3004 2.6214 3]; 
原 > 程 组 Ah 价 \/ 标 准 形 式 XP=[0:0.1:1.0]; 
旋 程 组 村 从 于 以 下 杜 准 形 z plot(xP,yP,"o',Xx,y(1,),r-) 


求解 : weliEeiitztsietslilieis 
solinit = bvpinit(linspace(0,1,10),[0 1]); 
eeetCL@]81 SLCLIGILTLIRESeliille 





的 方程 组 : [SteltATEW Te|ESIellniTeNLITTIT3Te1ESIelIniTeial[erectireiN[eialt rs 
[elelleexkeiil 
2 
[机 function dydx = ODEfun(X,y) 
dydx = [Y42); (1+X^2) y(1)+1]; 
人 
边界 条 件 为 : function bc -= BCfun(yayb) 


bc = [ya(1)-1;yb(1)-3]; 


边 值 问题 的 求解 








welteiitsietsliiie7 
求解 : c=1 |; 
[ielllileEaexheliilittiSiertes(e 区 [及 让 站 
[eaeeie(CL@1B1 SiC 1G31ELRESiel| alla 区 e 
X = [0:0.1:4]; 
y = deval(Sol,X); 
plot(X,y(1,:),b- ,sol.X,Sol.y(1,),，ro)) 
legend( 解 曲线 ,初始 网 格 点 解 ") 
2%6 


function dydx = ODEfun(X,y,C) 


他 





function bc = BCfun(ya,yb) 
bc = [ya(1);yb(1)+2]; 





边 值 问题 的 求解 辐 








Teiel 攻 GTeieltiilielsl 


要 ea 
求解 : Zz+(4-2qgcos20z=0 solinit = bvpinit(linspace(0,pi,10),[1 人 ,Imb); 


jelkaeMeiattsliciSeliD 
| [eee ie(CL@1B1 LEIRCD siG310LRESiel| il 下 ee 
lambda=sol.parameters 
边界 条 件 : x = [0:pi/60:pi]; 


y = deval(Sol,X); 
Z0)=1Lz(0)=0z(D=0 end 名 量 级 和 多 网 入 上 用 

legend( 解 曲 线 , 初始 网 格 点 解 ) 

function dydx = ODEfun(x,y,Imb) 

e 忆 5 

0 = [y(2); -(Imb-2 9q cos(2 X))Y(1)]; 


人 bc = BCfun(ya,yb,Imb) 
= [ya(1)-1;ya(2);yb(2)]; 


本 例 中 ， 微 分 方程 与 参数 入 的 数值 有 关 。 一 般 而 言 ， 对 于 任意 的 入 值 ， 该 问题 
无 解 ， 但 对 于 特殊 的 入 值 ( 特 征 值 ) ， 它 存在 一 个 解 ， 这 也 称 为 微分 方程 的 特 
征 值 问题 。 对 题 ， 可 在 bvpinit 中 提供 参数 的 铺 测 值 ， 然 后 重复 来 解 BVP 
得 到 所 需 的 参数 ， 返 回 参数 为 sol.parameters 


边 值 问题 的 求解 喇 


LIwlelReittseisiiets 
SelleEaeellaltdSieztws( 央 民 3) 太 Ci 
jelleliSEaeersisiittsitciSNei LEel Lo 

[ie 有 eeie((CL@1B1 ai 21G310LRESiel| aliEejeiileliiS3 
X = SO|.X; 











function dydx = ODEfun(x,y) 


dydx = [0.5"y(1)(Yy(3) - Yy(T))MY(2) 
边界 条 件 : -0.5*(y(3) - y(1)) 
(0.9 - 1000"(Yy(3) - y(95)) - 0.5"y(3)Y(G3) - YCTD)))Y(4) 


0.5 (Y(3) - y()) 


100*(y(3) - y(5)) ]; 
初始 猜测 解 为 站 


function bc = BCfun(ya,yb) 
bc = [yal(1) -1 
yatEj | 
ya(3) - 1 
) 
) 


function v = ex1init(X) 
=[1; 1; -4.5*X^A2+8.91*Xx+1; -10; -4.5"X^A2+9"Xx+0.91 ]; 








边 值 问题 的 求解 


Ce 








Infinity = 6; 

[el lllleEalexheliilittlarSierztesx(eLilalwAs) 居 LN 
jelleliSEaeeSisiitsiETESINeLDD 

[ie 有 eere((CI=》ksiele[ = 页 (CI=》《slele 丙 ie] 中 ojelilelitS 

eta = SO|.X; 

f = Sol.y; 

fprintf( An ); 

fprintf('Cebeci & Keller report f”"(0) = 0.92768、n') 

fprintf( Value computed here is f"(0) = %7.5fAn',f(3,1)) 

clf reset 

plot(eta,f(2,)); 

axis([O infinity 0 1.4]); 

title(Falkner-Skan equation, positive wall shear, \beta = 0.5.) 
xlabel( eta ),ylabel(dfdeta ),shg 

26 -7 
LIeteliEelielsitc 本 《sielel= (al 已 朋 

beta = 0.5; 

dfdeta = [f(2); f(3); -f(1)f(3) - beta*(1 - f(2)^2) ]; 


function res = ex5bc(fO,finf) 
res = [f0(1); f0(2); finf(2) - 1]; 





常 微分 方程 边 值 问题 的 应 用 
已 知 在 球形 氧化 铝 众 化 剂 进行 环 己 烷 脱 所 反应， 催化剂 直 径 为 


5mm， 操 作 温 度 为 700K， 在 此 温度 下 ，k = 4S-1，D=0.05cm2/s.。 
计算 众 化 剂 颗 粒 内 环 己 烷 的 浓度 分 布 及 俯 化 剂 的 有 效 因 子 。 


数字 模型 : 


0 











LIwieliEGittsieisiile 2 
global fai Cs 
Cs =1;rp =5e-3/2;k = 4; D = 0.05e-4; fai = rp*sgqrt(WD) % fai: Thiele 模 数 
a=0b=1; 
solinit = bvpinit(linspace(a,b,100),[0 0]); 
sol = bvp4c(@ODESs,@BCfun,solinit); 
plot(sol.x,sol.y(1,),b-) 
[el =] [人 BITITL=TIL [eictelIL 动 Eleis[ 人 BILLSieiil[ss Geiliwsttileliiy 
11 = quadl(Q@IntFunc1,a,b, 由 ,Usol.x,sol.y(1 ,2)) 
[ee (CLII 有) 
eta = |1/|12 
fprintf(\t 众 化 剂 的 有 效 因 了 于 为 : n = %.4fn',etal) 
526 人 
LielEeNe[E GD (RM 
global fai Cs 
dy1dr = y(2); 
让 丰 三 三 

dy2dr = 0; 
else 

dy2dr = fai^2*y(1)/Cs - 2/ry(2); 
ie 


function bc = BCfun(ya,yb) 
bc = [yb(1)-1; ya(2)]; 


function f= IntFunc1(rriyi) 
f= spline(riyi,m frA2; 


function f= IntFunc2(n) 
global Cs 
f = Cs” fA2; 





化 工 数 学 模型 中 的 微分 方程 Ge 


1. 集中 参数 模型 和 分 布 参 数 模型 
集中 大 数 模型 是 # 因 变量 不 随 空间 坐标 变化 ， 模 型 中 自 变量 仅 为 时 
间 ; 分 布 参 数 模型 则 指 因 变 量 不 仅 与 时 间 有 关 ， 而 且 与 空间 有 关 。 集 
参数 模型 为 党 和 ws 方程 。 

2. 人 常见 的 常 微分 方程 边 值 问题 模型 为 涉及 传 热 和 
扩散 的 问题 

3 在 建立 微分 方程 模型 时 ， 请 注意 初始 和 边界 条 件 ， 这 样 才 是 
一 个 完整 的 定 解 问题 

4 化工 模 型 中 的 微分 方程 种 类 很 多 ， 但 其 求解 过 程 往往 并 不 复 


习 。 


第 六 讲 
偏 微 分 方程 数值 解 


化 工学 院 软 件 应 用 教科 组 
2006-10 


本 章 知 识 要 点 





”MATLAB 求 解 侦 微 分 方程 图 数 一 pdepe 的 使 用 
方法 

。 [Matlab 的 PDE Toolbox 图 形 用 户 界面 CGUI) 
求解 侦 微 分 方程 

。PDE Toolbox 命 令 行 编程 求解 侦 微 分 方程 


偏 微分 方程 在 化 工 模型 中 的 应 用 





* 传 热 、 传 质 过 程 主要 是 抛物 型 方程 
椭圆 型 方程 ( 稳 态 方程 ) 
“流体 流动 方程 多 数 为 双 曲 型 方程 





偏 微分 方程 分 类 < 


线性 偏 微分 方程 : 
方程 经 过 有 理化 并 消去 分 式 后 ， 若 方程 中 没有 未 知 函 数 
及 其 偏 导数 的 乘积 或 圭 等 非 线 性 项 ， 那 就 称 之 为 线性 偏 
微分 方程 (组 ) 。 如 果 有 则 称 为 非 线性 偏 微 分 方程 

拟 线 性 偏 微分 方程 : 
如 果 仅 对 未 知 函 数 的 所 有 最 高 阶 导 数 来 说 是 线性 的 ， 则 
称 之 为 拟 线性 偏 微分 方程 (组 ) 。 

齐 次 与 非 齐 次 方程 : 
在 线性 方程 中 ， 不 含有 未 知 函 数 及 其 偏 导 数 的 项 称 为 自 
由 项 ， 自 由 项 为 零 的 方程 称 为 齐 次 方程 ， 否 则 称 为 非 齐 


次 万 程 。 


二 阶 偏 微分 方程 的 分 类 








当 式 中 A，B，@C 为 常数 时 ， 即 为 拟 线性 偏 微分 方程 ， 可 以 分 为 三 类 : 


当 YU ij ， 方 程 称 为 圆 型 (elliptic ) 偏 微分 方程 ; 
当 EEE CE 人 ij ， 方 程 称 为 抛物 型 (parabolic ) 偏 微分 方程 ; 
当 PR YE ij ， 方 程 称 为 双 曲 型 (hyperbolic ) 偏 微分 方程 。 





三 个 最 基本 最 典型 的 偏 微分 方程 


双 曲 型 方程 (波动 ) 方程 





椭圆 形 方程 ( 拉 普 拉 斯 方程 ) 





抛物 型 方程 (热传导 方程 ) 








偏 微分 方程 的 定 解 条 件 
初始 条 件 : 攻 有 


对 应 初始 条 件 的 定 解 问题 称 为 柯 西 (Cauchy ) 问题 








对 应 的 定 解 问题 称 为 Dirichlet 问 题 
对 应 的 定 解 问题 称 为 Neumann 问 题 


人 only 对 应 的 定 解 问题 称 为 Robin 问 题 


既 有 初始 条 件 又 有 边界 条 件 的 定 解 问题 称 为 混合 问题 








偏 微分 方程 的 数值 解法 








有 限 差 分 法 : 
微分 方程 的 导数 用 有 限 差 分 近似 代替 ， 最 后 得 到 代数 方程 组 
求解 偏 微分 方程 最 基本 的 方法 

正 交 配置 法 : 
采用 正 交 多 项 式 为 试 函 数 ， 构 造 其 系数 使 其 近似 于 微分 方程 近似 解 
求解 边 值 问 题 的 有 效 方法 

线 上 法 : 
将 一 个 自 变 量 当 成 连续 变量 ， 而 对 其 自 变 量 用 有 限 差 分 法 或 正 交 配 
置 法 进行 离散 ， 从 而 把 偏 微分 方程 转变 为 常 微分 方程 组 

有 限 元 法 : 

变 分 原理 和 分 片 播 值 的 结合 体 ， 将 求解 区 域 分 成 有 限 个 单元 ， 然 后 
采用 分 片 播 值 函 数 逼 近 解 ， 从 而 得 到 代数 方程 组 。 


线 上 法 (MOL ) 求解 偏 微分 方程 示 爷 


00, “0,) = 0 


。 ”在 Xx 方 向 ， 将 求解 区 域 5 等 份 ， 则 下 的 值 可 用 中 间 四 个 节点 
处 的 值 表示 

。 下 对 x 的 2 阶 导 数 以 差分 代替 ， 即 医生 且 

。 ， 原 方 程 求解 转化 为 以 下 常 微 分 方程 组 的 初 值 问题 











Matlab 求 解 偏 微分 方程 方法 


pdepe 马 数 
人 求解 一 维 动态 模型 ， 即 一 维 抛物 型 和 椭圆 型 偏 微分 方程 的 初 值 - 边 值 
问题 


他 采用 线 上 法 求解 ， 利 用 正 交 配置 进行 空间 离散 


pdetool 图 形 用 户 界 面 
他 求解 二 维 线性 和 非 线性 偏 微分 方程 的 PDE 工 具 箱 
他 有 限 元 法 
他 可 用 于 求解 具有 标准 形式 的 偏 微分 方程 


pdepe 求 解 偏 微分 方程 的 形式 










偶 微 分 方程 组 


m = 0，1，2 分 别 代表 平板 ， 圆 柱 和 球形 








pdepe 马 数 的 贡 用 





调用 格 亏 : 

ie elelsels0leielse 同 le 本 Ts Ge 

Sie 呈 eielseis0leielse 同 wiiwilT GeceleliS 

Sie 吴 eielaes0l 本 else 汪 wliewilis Ge 本 eleii 二 


输出 参数 : 

输出 变量 Sol 是 三 维 数组 

急 Sol(:,:;)=ui 是 解 向 量 的 第 i 个 分 量 

急 Sol(,k,)=Uu(,k) 是 解 向 量 Ui 在 (xb=(tspan(),xmesh(k)) 处 的 数值 解 

急 Sol(j,:, 站 是 解 向 量 第 i 个 分 向 量 uj 在 时 间 tspan() 和 网 格 点 Xmesh(:) 处 的 数值 解 





pdepe 函 数 的 输入 参考 好 





入 参数 : 


m 表 示 定 解 问题 的 对 称 性 ，m=0, 1, 2 分 别 代表 平板 、 圆 柱 和 球形 。 

pdefun 是 描述 PDE 问 题 的 函数 ， 输 出 变量 为 标准 形式 中 的 c,f 和 s， 其 格式 为 [c,f， 
s]=pdefunc(x,tiududx)。 输 入 变量 X 和 t 为 标量 ，U 和 dudx 是 向 量 。U 和 dudx 分 别 是 
问题 的 解 U(XxU 和 它 对 X 的 偏 性 数 的 近似 。c,f, s 是 列 向 量 。 

icfun 描 述 定 解 问题 初始 条 件 的 函数 ， 其 格式 为 U=icfun(X)。icfun 计 算 和 返回 解 的 初 
始 值 。 

bcfun 是 描述 定 解 问题 边界 条 件 的 函数 ， 格 式 为 [pl, ql, pr, qrl=bcfun(xl， ul, xr，ur, it)。 
其 中 ul 是 在 左边 界 Xxl=a 处 的 近似 解 ，uUr 是 在 右边 界 xrf=b 处 的 近似 解 。pl 和 ql 是 Pp 和 q 
在 X| 处 的 列 向 量 值 ， 同 样 pr 和 qr 是 p 和 qd 在 xf 处 的 列 向 量 值 

xmesh 是 空间 变量 x 的 网 点 向 量 ，pdepe 不 会 自动 选择 。xmesh=[x0, x1, ..…… xn]， 
满足 X0<X1<...<Xxn， 且 xmesh 的 长 度 必 须 大 于 3，xmesh(1)=a 和 xmesh(end)=b。 
pdepe 的 求解 效率 与 Xmesh 的 选择 好 坏 关 系 很 大 。 通 常 在 梯度 较 大 处 应 加 密 网 格 。 
tspan 是 时 间 变 量 t 的 网 点 向 量 ， 是 用 户 在 该 处 想得到 数值 解 的 时 间 向 量 。tspan=[t0， 
t1,，..., tn]， 满 足 t0<t1<...<t， 且 tspan 的 长 度 必 须 大 于 3，tspan(1)=t0 和 
tspan(end)=tf。 求 解 效 率 和 tspan 的 牙 密 关 系 不 大 。 





pdeval 允 数 的 使 用 


。 用 于 计算 在 区 间 [a,pb] 内 的 函数 值 和 对 X 的 偏 必 数值 

。 调 用 格式 : [xoutduoutdx]=pdeval(m, xmesh, ui, xouf) 

*Uji=sol(G ,iD 是 问题 解 的 第 ji 个 分 量 在 时 间 tt 和 网 点 xmesh 处 的 值 。 
。XOUt 有 是 [X0,xn]=[ab] 内 的 数 构成 的 向 量 

。uout 和 duoutdx 分 别 是 芭 本 和 ER 玖 在 xout 处 的 马 数 值 


一 维 热传导 方程 求解 


一 维 热传导 方程 





边界 条 件 : 


4 条 什 。 FEOREITC 








一 维 热传导 方程 求解 











原始 方程 初始 条 件 界 条 件 
县 上 县 





ER mm 
上 





一 维 热传导 方程 求解 


TelelEGittcielsiiiie 

m = 0; 

Xx = linspace(0,1,20); 

t = linspace(0,2,5); 

je] 省 ele[sjelsi( 区 Cole 人 ieie[= 丰 (Cole 和 TeCCJeie[ > 和 ee 入 和 和 
% Extract the first solution component as U.， 

U = Sol( 1); 

esieliemETeitcTiEE2Ee[eieleEWLcN 人 esiLULeN 人 :ESielieliE 
Surf(X,t,U) 
ENLITLTLSITe1IESeLUieIEeellileigitsieRtE20ILs ieieiiil 测 ) 
xlabel('Distance X") 

ylabel(Time tt) 

人 ell yiileliiielrelil ez1EELLieaei=| Ilattiliie 

figure 

plot(x,U(end,:)),title(CSolution att = 2) 

Xlabel(' Distance x ),ylabel(U(X,2)) 








function [c,fs] = pdex1pde(x,tiuDUDXI) 
cC=pi^2f=DuDxi:s = 0|; 

2p0 人 
function U0 = pdex1ic(X) 

U0O = Sin(piX); 


function [pl,ql,prqn = pdex1bc(xlul,xrurt) 
pl = Ul;ql= Opr= pi*exp(-b);qr = 1; 


偏 微分 方程 求解 





定 解 问题 标准 形式 








偏 微分 方程 来 解 





function Cha6demo2 

m = 0; 

X = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; 
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2]; 

sol = pdepe(m,@pdex4pde,@pdex4ic,@Opdex4bc,Xx,t); 

u1 = Sol(:,:,1);uU2 = Sol(:,:;,2); 

figure 

Surf(Xx,t,u1),title(Cu1(Xx,t)),Xxlabel(Distance Xx"),ylabel(Time t) 
figure 

Surf(X,t,uU2) 

title(u2(Xx,t)),xlabel(Distance Xx' ),ylabel(Time t) 


function [c,f,s] = pdex4pde(x,tu,DuDx) 
c= [1; 1]; 

f = [0.024; 0.17] ” DuDxX; 

y = U(T) - u(2); 

F = expb(5.73*y)-exp(-11.47*Vy); 

s = [-F; F]; 


function [pl,ql,prqr] = pdex4bc(xlulxrurt) 
pl = [0; uUI(2)]; 

ql = [1; 0]; 

pr= [ur(T)-1; 0]; 

qr = [0; 1]; 








急 在 命令 窗口 键入 pdetool， 则 进入 如 右 图 所 示 的 图 形 用 户 界 
面 (GUI ) 


四 Lel zl sx 








Rile Bait Dptions Dear Bounaary TIE Mesh Solve 了 Lot 量 ndow elp 
| | | | | | | 
| 司 | 全 | 二 | | |Poe| 全 | 念 | = 4 CA | ceneicscar 加 | 国 sssssss 


Set formula 






























































GUI 能 求解 的 各 类 偏 微分 方程 形式 





椭圆 天 方 程 : 
抛物 型 方程 : 
双 曲 瑾 方程 : 
特征 值 方程 : 
非 线性 顶 圆 型 方程 : 


偶 微 分 方程 组 : 





GUI 求解 偏 微分 方程 的 定 解 条 件 


en 


”在 非 线 性 条 件 下 hpq,g 可 以 是 U 的 邓 数 ; 对 于 抛物 
型 和 双 曲 型 方程 而 言 ， 可 以 是 时 间 t 的 函数 


对 于 二 维 体系 : 
Dirlchlet: 





IN[s1WTE1TAE 


混合 边界 条 件 : 





GUI 求解 偏 微分 方程 的 步骤 


含 采 用 GUI 求解 偏 微分 方程 ， 其 步骤 如 下 : 


1 
.设置 边界 条 件 


co ~1oDwO 上 mm 


画 求解 区 域 


设置 方程 


. 网 格 剖 分 


解 方 各 
图 形 结果 


.输出 网 格 
.输出 解 的 数值 





GUI 界面 的 菜单 





. File 

Edit 基本 操作 
本 @jeiielii 
Draw 
ielllileic1 
PDFE 

. Mesh 
.Solve 

9. Pilot 
AIele 
本 je 


oo ~lIODO 上 NOmND 一 


品 | 田 | 全 | 二 | 闻 appPog 人 人 一 人生 





Enters the boundary mode. 


Opens the PDE Specification qialog box. 


Initializes the triangcular mesh. 


Refines the triangular mesh . 


Solves the PDFE. 


3-D solution opens the Plot Selection dialog box. 











霄 了 -EDE Ioolboz = [Untitledj] 同 别 咱 霹 可 引 区 | | 


了 il 是 File 了 ait Options Draw ”Boundary FDIE Mesh Solve 粳 indow ”Help - 
人 了 | X， -0.6115 EREHETO2 


[上品 品 下 le[z|el 天 AI 全 | = 条 反 















0.15 


EC 
EN 


0.05 





-1.5 -1 0.5 D 0.5 1 1.5 
Info: Pushto initialize mesh . 
轴 nanonoon | 二 








GUI 求解 抛物 型 方程 示例 


热传导 方程 : 金属 板 导 热 问 题 
一 个 带 有 适 形 孔 的 金属 板 ， 板 的 
左边 保持 在 100C， 右 边 热 量 从 板 
向 环境 空气 定常 流动 ， 其 它 边 及 
内 孔 边 界 保 持 绝热 。 初 始 时 板 的 
温度 为 OYC 。 

金属 板 的 外 边界 坐标 为 (-0.5,- 
0.8)，(0.5, -0.8), (0.5,0.8)， (- 
0.5,0.8)， 内 边界 坐标 为 (-0.05，- 
0.4), (-0.05, -0.4), (-0.05, -0.4),(- 
0.05, -0.4) 





Time=5 Color uU Height: U 





Au=00yeo={ey1<x<L-I<y<]| 














]C 0.5 
&U=arctali co 习 已 三 才 


-0.5 


宇 -3sinboge 207 二 0 
Cr 
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收 是 





1 日 广 厂 EEC 


ee 


全 放 分 为 多 个 子 区 域 人 ， 在 每 个 单元 里 以 
种 简单 的 函数 来 代替 未 知 函 数 
多 如 果 0 是 方程 的 解 以 任意 的 试 函 数 并 积分 : 





售 上 式 分 步 积 分 ， 并 代入 边界 条 件 ， 则 原 求解 问题 变换 为 如 
下 问题 : 求解 U 满 足下 式 : 


(er 





有 限 元 法 的 基本 思想 


人 
| 


ER 风度 和 件 


0 
国 7 
加 了 


抛物 方程 ae :KU - 王 ED 














Boundary 
Condition 
matrix 


Geometry 
Description 


matrix 
decsg 
Decomposed 
Geometry 
matrix 
initmesh 
Boundary Mesh 
M-file data 
assempde 
Solution 
data 


pdeplot 


Geometry 
M-file 


refinemesh 


Coefficient 
matrix 





Coefficient 
M-file 





PDE Toolbox 马 数 :用 户 界 面 算 法 


目的 


eleliillie 瑟 李 圆 


转化 为 1.0 版 本 的 M 文 件 


ele[eey 国 多 边 形 
进入 工具 箱 的 用 户 界面 











PDE Toolbox 马 数 : 几何 算法 


分 解 结构 矩阵 成 最 小 区 域 
创建 初始 网 格 

十 

加 





在 矩形 区 域 上 创建 规则 网 格 
可 三 角形 Id 属 


写 边 界 条 件 说 明文 件 
写 几 何 说 明文 件 








PDE Toolbox 函 数 : PDE 数 值 计 算 函 数 


组 装 积分 区 域 贡献 ，K,M,F 
组 装 刚度 矩阵 和 方程 右边 的 函数 
蛙 双 曲 型 PDE 问 题 


和 而 PE 
fiPDEN 


非 线性 PDE 问 是 
求 矩形 网 格 上 泊 松 方程 的 快速 





PDE Toolbox 允 数 : 


绘图 函数 用 户 定义 算法 
eleleelil Nt esieiellile 由衣 这 值 线 国 的 
速 命 令 速 命 令 
IPDE 几 何 图 全 PDE 儿 何 图 





绘制 PDE 三 角形 
网 格 册 


pdeplot | 一 般 PDE 工 具 箱 
会 锋 国 数 


会 制 表面 图 的 速 





PDE Toolbox 允 数 : 工具 算法 


用 相对 误差 判别 准则 选 子 区 域 集中 边缘 的 索引 
择 最 低劣 的 三 角形 了 区 工人 中 点 的 家 
相对 最 差 信 远 择 低 和 
本 $ 构 力学 张 生 了 




















三 角形 几何 数据 
有 邻 的 三 三 角形 几何 数据 


角形 的 索引 测量 三 角形 网 格 质量 


从 单元 节点 数据 到 三 角 | sptarm 解 一 般 稀 玻 矩阵 特征 值 
形 单元 中 点 的 插值 问题 


从 三 角形 中 点 数据 到 节 | tri2grid 由 PDE 三 角形 网 格 转换 
点 数据 的 插值 成 矩形 网 格 











命令 行 求解 PDE 的 一 般 步 又 好 


1. 问题 定义 及 参数 初始 化 - pdetool wbound, wgeom 
定义 求解 区 域 (几何 条 件 ) 、 边 界 条 件 和 方程 系数 及 参 
数 初 始 化 
2. 网 格 化 (三 角形 网 格 划 分 、 网 格 细 化 及 微调 ) -initmesh， 
refinemesh, lgglemesh 
为 便于 观察 网 格 ， 可 以 采用 pdemesh(p,e 让 绘制 由 网 格 
数据 p,eit 制 定 的 网 格 
汪汪 尘 和 BT-11jeret- eic|eieleaaieisieel[eeie[siiielilii 
ee 本 ee[sewelili 


请 通过 Help 命 令 查看 以 上 命令 的 调用 格式 


了 


命令 行 求 解 椭圆 方程 示例 





求解 单位 圆 上 的 泊 松 方程 : 








gs 求解 问题 的 几何 和 边界 条 件 已 分 别 由 circleg， ECTSDTL 定 义 


[p,e,l=initmeshCcircleg , Hmax ,1); 一 初始 化 网 格 
[p,etl=refinemesh(circleg ,p,e,t); 

[p,e,tj=refinemesh(circleg ,pe,b; 网 格 加 密 
[p,e,tl=refinemesh(circleg ,p,e,t); 


WEISSieeiGeliels ee 有 机 本 0 区 中 


Ai 一 一 一 一 绘制 网 格 图 


pdesurf(p,tu) ”绘制 解 曲 面 图 








命令 行 与 GUI 联合 求解 偏 微分 方程 示例 


求解 单位 圆 上 的 泊 松 方程 : 


1 . 
2. 


3， 








利用 GUI 定义 几何 模型 和 边界 条 件 ; 

在 GUI 的 [Boundary] 菜 单 中 ， 选 择 [Export Decomposed Geometry， 
Boundary Cond's...]， 将 几何 模型 和 边界 条 件 输出 

人 eellliie[(e 区 elel[sisiellle 翅 有 2 elselilife 同 eleltSseliie1 
将 几何 模型 和 边界 条 件 存 入 文件 poissong.m 和 poissonb.m 


. 采用 以 下 语句 求解 方程 并 输出 结果 


EeelSiieliie 刘 

6 二 elel[siielile 冰 

c=1;a=0:f=1; 
[p,e,tj=initmesh(g)， 
[petj=refinemesh(g,p,e,t); 
U=assempde(b,p,etic,a,f); 
pdesurf(p,t,U) 
(welleliiitciel(eele]j 


等 温 一 级 反应 的 圆柱 状 催化 剂 颗粒 中 浓度 分 步 加 肛 


号 了 


function CylindCat1 
eeLGi 

b = 'CylindCat1b ; 

C = 'X;a=0;f= -4 XU 

[pe = initmesh(g) 

[p,eil = refinemesh(g,p,e,b; 

p = jigglemesh(p,e,tb); 

[ures] = pdenonlin(b,p,e,tic,a,f) 
ele[ssigLiIge 有 机 看 碎 eelleller1l 
xlabel(Xx'),ylabel(y'y),zlabel(OuU ) 











利用 Matlab 求 解 偏 微分 方程 的 策略 





对 于 一 维 动态 模型 ， 优 先 使 用 Matlab 函 数 pdepe 
2. 对 于 二 维稳 态 和 动态 模型 
。 含 有 两 个 方程 以 下 的 标准 PDE 问 题 ， 优 先 采 用 pdetool 
。 含 有 三 个 方程 以 上 的 标准 PDE 问 题 ， 选 择 PDE 工 具 箱 
3. 当 以 上 方法 均 无 法 解决 问题 ， 考 虑 使 用 正 交 配置 法 ， 有 
限 差 分 等 方法 。 对 于 许多 化 工 偏 微分 方程 问题 而 言 ， 
MOL 法 也 是 不 错 的 方法 


第 七 讲 
应 用 数理 统计 和 了 最 优化 方法 


化 工学 院 软 件 应 用 教科 组 
2006-11 


化 工 数 学 模型 建立 的 一 般 步 骤 





最 优化 方法 ， 
ED 
2 数理 统计 | 





本 章 知 识 要 点 





。 Matlab 统 计 工 具 箱 函 数 
”Matlab 最 优化 工具 箱 函 数 

。 Matlab 遗 传 算法 与 直接 搜索 工具 箱 
。 最 小 二 乘法 线性 、 非 线性 曲线 拟 合 





Statistics Toolbox 的 内 容 分 类 


0 人 elelliiy 人 Bidleigiiel 

。 描 述 统计 (Descriptive Statistics ) 

。 统 计 可 视 化 (Statistical Visualization ) 

“假设 检 验 (Hypothesis Tests ) 

。 线 性 模型 (Linear Models ) 

。 非 线性 回归 (Nonlinear Regression Models ) 
。 多 元 分 析 (Multivariate Statistics ) 

。 统 计 的 过 程控 制 (Statistical Process Control ) 
。 实 验 设 计 (Design of Experiments ) 

。 隐 马尔 科 夫 模型 (Hidden Markov Models ) 





需 机 变量 与 概 府 分 布 Ce 
as 
9 变量 


。 概 率 分 布 : 人 人 弘 属 于 一 给 定 值 
集 的 概率 所 确定 的 函 


-分 布 了 数 : Da 
0) = FOOD) 1 的 
_ 分 布 函数 的 逆 画 数 也 各 为 分 位 数 : 上 于 末 








Statistics Toolbox 中 的 概率 分 布 种 类 


*Maitlab 提 供 约 30 种 概率 分 布 ， 主 要 的 如 下 : 
2 人 elliteliili1 区 elisiidlelyiieli 
1 条 ielalif| eltsiidleiyiiei 才 
2 并 作 elejiseliEeiisildlelii[ell 
2 信人 11 区 :eliTYeltsildeiliiell 
一 t 分 布 全 :ie[ilwge[siidleniiie 人 
一 F 分 布 生 e|iSildleiliiel 
eleli[c1iiii[EelisidleIniiei 二 





概 座 分 布 的 功能 函数 


。 对 于 每 种 概率 分 布 ，Matlab 提 供 了 如 下 6 种 功能 函数 
ee] 
(edijeiel 本 wiiWeiel 辣 etcieel 明 elilejeieliaieeljeelR 间 
ee 
(elwel 本 il 放 we 天 ectweliellieeellteel 本 we 
一 分 位 数 函 数 ; “inv 
一 分 布 统计 函数 ; “stat 
一 分 布 拟 合 函数 册 宙 | 
一 随机 数 生 成 函数 ; “rnd 





正 态 分 布 的 计算 


系 积 分布 函 数 normcdf: 
例 : 标准 正 态 分 布 在 区 间 (-1,1) 内 的 取 值 概率 

p=normcdf([-1 1],0,1); P(2)-p(1) 

系 积 分 布 函 数 normpdf: 

例 : 标准 正和 态 分 布 ,X=1 的 概率 客 度 

p=normpdf(1;0,1) 

ed 

例 : 标准 正 态 分 布 的 包含 95% 取 值 的 区 间 

Xx=norminv([0.025 0.9751,0,1) 


可 通过 disttool 打 开 分 布 函 数 图 形 界 面 ， 观 察 计 算 结 果 


描述 ,性 统计 


。 位 置 特 征 
一 算术 平均 ，mean 
一 中 位 数 ，median 
一 切 尾 平均 ，trimmean 
一 几何 平均 ，geomean 
-调和 和 平均 ，harmmean 


一 众 考 ) Jelel: 


。 变 措 特 征 


-- 极 差 ，range 
-平均 绝对 偏 关 ，mad 
1 
re 
Te 








假设 检验 
。 假 设 检验 是 一 种 为 了 确定 一 个 有 关 总 体 分 布 的 假 
设 应 该 子 以 绝 还 是 接受 的 程序 。 
。Statistics Toolbox 假 设 检验 函数 包括 以 下 两 类 : 


“统计 检验 "分 布 检验 
一 dwtest 一 Chi2gof 
一 ranksum piesi 
一 runstest 0 
-signrank -lillietest 
eliltsisii 
一 ttest 
一 ttest2 
一 Vartest 
一 Vartest2 


一 ztest 





。 多 元 线性 回归 函数 : regress 
b = regress(y, 久 ) 
[b,bint] = regress(y, 义 ) 
[b,bint,r] = regress(y, 久 ) 
[b,bint,rrint] = regress(y, 义 ) 
[b,bint,rrint,stats] = regress(y, 久 ) 
[b,bint,r,rint,stats] = regress(y,X,alphal) 
其 中 : 
b， 回 归 系 数 ; bint， 回 归 系 数 的 置信 区 间 
r， 残 差 向 量 ; rint， 残 差 向 量 的 置信 区 间 
stats 包 含 四 个 值 : 相关 系数 、F 统 计量 值 相 应 于 F 统 计量 值 的 概率 、 误 差 方 差 
估计 
alpha， 置 信 水 平 





regress 的 使 用 


已 知 某 观测 量 y 的 值 与 X1,X2， 
= 个 恋 量 有 关 ， 如 左 表 上 

XO3 三 1 变量 有 》 如 左 表 所 
区 
示 ， 采 用 线性 同 轨 裤 合 邓 数 
y=[4.3302 3.6485 4.4830 5.5468 5.4970 3.1125.… 
四 5 2 5.0060 5.2701 5.3772]'; 

人 | 一 ! 

x2-[189143201217582316184]; 
3.8759 |6 ”15 |1s9 | xs-[50404643644064393755604950] 


Xx0=ones(1,13); 









FE 
5.0060 四 靖 6 6 60 
5 
5 





多 项 式 回归 


。 多 项 式 回归 的 主要 函数 
1.[a,Sl=polyfit(X,y,m) 
2.[y,DeltaAj=polyval(a,X,S) 
3.[y,DeltaAj=polyconf(a,x,S,alphal) 
。 如 果 输 入 数据 的 误差 相互 独立 ， 且 有 常数 方差 的 正 态 
分 布 ， 则 yt+ DeltaA 包 含 100(1-alpha)% 的 置信 区 间 
人 
。 给 出 用 n 阶 多 项 式 拟 合 X,y 数 据 ， 并 用 描绘 100(1- 
alphaje6 置 信 区 间 


除 线性 回归 和 多 项 式 回归 外 ， 统 计 工 具 箱 中 线性 回归 的 分 类 中 还 
包括 二 次 响应 面 模型 、 广 义 线性 回归 、 逐 步 回 归 





线性 模型 - 方差 分 析 


。 方 差分 析 (Analysis of Variance,ANOVA) 是 研究 一 
些 因素 〈 自 变量 ) 对 茶 个 指标 〈 因 变量 ) 的 相关 
关系 ， 研 究 因素 对 指标 的 影响 是 否 显著 


。 单 因素 实验 方差 分 析 函 数 : anoval 
。 双 因素 实验 方差 分 析 函 数 : anova2 
。 多 因素 实验 方差 分 析 函 数 : anovan 


多 元 统计 分 析 


主 成 分 分 析 : 考虑 多 个 数值 变量 间 相 关 性 的 一 种 
多 元 统计 方法 

eeeliileNNeiecUwellNeiecile ecLtcIsil 

聚 类 分 析 : 多 元 统计 中 研究 如 何 " 物 以 类 聚 " 的 统 
计 方 法 

-clusterdata，pdlist，sdquareform，lInkage 

判别 分 析 : 研究 如 何 根据 已 知 归 属 类 型 的 一 批 数 
据 来 推断 未 知 归属 类 型 的 数据 中 每 个 个 体 所 属 类 
型 的 一 种 多 元 统计 方法 

ENNULEULEI 

多 元 方差 分 析 

LE1Le 本 UETeUeiISic 1 





实验 设计 


。 完 全 析 因 设计 (Full Factorial Desing ) 
一 fullfact，ffe2n 

。 不 完全 析 因 设计 (Fractional Factorial Design ) 
一 [fracfact 

。 响 应 面 设 计 (Response Surface Design ) 
lele[s :jie wwetssiell 

9 吧 @jeiiliit| 四 BITi[eli 才 


一 manova，manovacluster 





最 优化 问题 的 一 般 形 式 

















目标 函数 线性 优化 
非 线 性 优化 
多 目标 优化 





Matlab 最 优化 函数 分 类 


求 极 小 值 问 题 
国 非 线性 最 小 化 CNonlinear minimization of functions ) 
图 多 目标 非 线 性 最 小 化 (Nonlinear minimization of multi- 
elej[ se ieieii 才 
团 最 小 二 乘 问题 
@ 线性 最 小 二 条 (Linear least squares of matrix problems) 
@ 马 数 的 最 小 二 乘 (Nonlinear least squares of functions) 
所 人 二 etieliEseliie) 
@ fzero,fsolveA、 


优化 函数 控制 





极 小 值 问 题 函 数 


整数 规划 
多 目标 达到 问 


有 边界 的 单 变 量 非 线性 最 小 化 问题 


有 约束 非 线性 最 小 化 问题 





二 元 整数 规 世 
可 题 





:放电 
无 约束 非 线性 最 小 化 
无 约束 非 线 性 最 小 化 
半 无 限 问 题 


二 雇 规 划 问 题 





函 数 fminbnd 和 单 变 量 最 优化 问题 








调用 格式 : 
x={fminbnd(fun,Xx1,X2) 
Xx=fminbnd(fun,x1,x2,options,….p1,p2) 
[x,fvall = fminbnd(.…) 
[x,fval,exitflag] = fminbnd(.…) 
[x,fval,exitflag,outputj = fminbnd(.…) 


输入 参数 : 

fun 目标 函数 ， 可 以 由 M 文 件 定 义 ， 或 由 匿名 函数 确定 
Xx1，X2 求解 区 间 [x1,x2] 

options 优化 参数 选项 


p1，p< 需要 传递 的 参数 


也 数 fminbnd 和 单 变量 最 优化 问题 








输出 参数 : 

X 最 优 解 

fval 目标 函 数 的 最 小 值 
iT 描述 退出 条 件 : 


Et 表示 目标 函数 收敛 于 解 X 处 ; 
exitflag = 0 表示 已 经 达到 函数 评价 或 和 迭代 的 最 大 次 数 ; 
exitflag<0 表示 目标 函数 不 收敛 


output 是 一 个 结构 体 ， 输 出 优化 过 程 的 各 种 信息 
elliein Lee 使 用 的 算法 ; 
output.func 马 数 评价 次 数 
output.iterations 和 迭代 次 数 


output.message 退出 信息 





甩 数 fminbnd 的 使 用 





在 区 间 [0 100] 内 求 取 单 变量 函数 的 最 小 值 





LIetelEGIiTXeL=liie 

Xx1 = 0; 

24 

[x,fval] = fminbnd(@2ObjFunc,x1,Xx2); 
fprintf( \nResults、n ) 


fprintf( Optimum solution: %fn ,X) 
fprintf(” Objective value: %f,fval) 
7 


function f = ObjFunc(X) 
f{ = XA3 + 3"XA2 - 9 X; 


马 数 linprog 和 线性 规划 











调用 格式 : 
人 erere[f 本 < 有 e) 
x=linprog(f,A,b,Aeq,beqd) 
人 elgere[( 到。 明 e 再 < ee 四 ee 加 le 及 le 及 48 黄 ejelt[eliis 
[x,fval] = linprog(.…) 
[x,fval,exitflag,output,lambda] = linprog(.…) 


其 中 lambda 是 一 个 结构 体 ， 人 钨 有 在 解 X 处 使 用 的 Lagrange 冬 
子 。 其 余 与 前 类 似 


函 数 linprog 的 使 用 








LewieliEEGTXeLsiile 灾 
f= [5 -6 -7]; 
A=[5 -6 10 
-1-5 3]: 
b = [20 -15]; 
re 
ee 
lb = zeros(3,1); 
NIHLETe 区 eglilell 本 Lele[1EE 二 人 [efere[h 丁 骨 e 克 Ace 四 el=ie 四 9) 


是 舍 











AX <_B min 5x1 - 6x2 - 7X3 
下 AeqX = Beq S.t. 5x1 - 6x2 + 10x3 <= 20 
LB<-X<-UB -X1 - 5Xx2 + 3Xx3 <=-15 
x{+Xx2+Xo=5 


Xx1,X2,X3 >=0 





马 数 fminunc，{fminsearch 和 无 约束 多 变量 问题 最 优化 


in700 本 


fminunc 和 fminsearch 分 别 采 用 的 是 拟 牛 顿 法 和 Nelder-Mead 单 纯 形 法 


< 一 


数学 模型 : 





调用 格 却 : 
x={fminuncl(fun,xuo,options...p1,p2) 
[x,fval,exltflag,output,grad,hesslan]| = fminunc(..) 
fminsearch 与 fminunc 调 用 格式 相同 。 





马 数 fmincon 和 多 变量 有 约束 最 优化 问题 (人 


一 


数学 模型 : 调用 格式 : 

x=fmincon(fun,x0,A,b) 

人 TeieliTiELI 生 《0 古本 e 再 AIe 四 ee 四 je 有 ie 用 [elileell 玉 eleillelatS 
[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(.…) 





台数 fmincon 的 使 用 








LileliEE@iiLTXeislilie7 
Xx0=[-1,1];a=[1,-1]j;b=1;a1=[1,1];b1=0' 
[x,fj=fmincon(@2objfun,x0,a,b,a1,b1,, 几 nonlconm) 


function f = objfun(XI) 
f = exp(X(1))(4*X(1)A2+2*X(2)^A2+4*X(1)X(2)+2*X(2)+1); 


ee [oa 区 e 汶 中 避 ellweiie4) 
c1=[1.5+Xx(1)*X(2)-x(1)-X(2);-x(1T)*Xx(2)-10]; 
C2=0; 





最 小 二 乘法 问题 


线性 最 小 二 乘 (Linear least squares of matrix problems) ) 
线性 约束 的 最 小 二 乘 
非 负 约束 的 最 小 二 乘 


函数 的 最 小 二 滋 (Nonlinear least squares of functions) ) 
最 小 二 乘法 曲线 拟 合 
证 最 小 二 灯 问 是 





lsqnonlin 可 以 来 解 问题 的 模型 如 下 : 





lsqcurvefit 可 求解 问题 的 模型 如 下 : 


See 人 si NeliielllilaisnzEe 








调用 格 式 : 
x= lsqnonlin(fun,x0,lb,ub,options) 
| 人 Lii[eldll 玫 Crsi[elg[c[ 区 =》>4liltcte 区 elglielg 明 [Leie[= 丙 [=Ueeleltc1]EE 二 和 ellellliiiI 人 本国 
Xx = |sqcurvefit(fun,x0,xdata,ydata,lb,ub,options) 
| 人 LateldllKsrsi[elg1I 区 siilLcTe 区 elglielgi 明 Lileie[ 再 [=Ueeleltc1]EE [ieeglasiir 全 办 


输入 参数 : 

fun 定义 F(X)， 注 意 函 数 返回 的 是 fun(X) 

X0 初 值 

[ee 设计 变量 的 下 边界 和 上 边界 ， 如 果 X(i) 下 无 界 ， 

则 设置 lb = -inf， 如 果 上 无 寞 ， 则 ub = inf。 

输出 参数 : 

X 拟 合 系数 

resnorm Xx 处 的 残 差 平 方 和 范 数 值 

residual X 处 的 残 差 值 


ELwelelEl Xx 处 的 Jacobian 拢 阵 


了 最 小 \ 二 乘法 非 线性 曲线 拟 合 


:sat: [aaETETETETSTTTE 


ouexTaicAxxus 国 瑟 3 











LewleiE Eee 
global xy 

XE= | 9， 

y=[15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6]; 
a0=[1,1]; 
[ERLsSLiLeldilRCssLelgtzi 刁 SielellAsiiit(C2xejejliili 区 10 及 的 人 
已 富 elileliliI(CielejlilPEsL0) 
Xx1=[0.2:0.6:8j]; 

y1=al(1) exp(a(2)X1) 
X2=[0:0.5:8]; 
y2=a2(1)exp(a2(2)X2); 
plot(Xx,y, bo ,xf1,y1,g.,Xx2;,y2,,m”) 


LIeleEelejiiyLiIK= 本 4) 
{=a(1) expPta(2)X); 


Lielllpee elejlilliT2k= 本 4) 
global x y 
for i=1:8 
zU)=a(1) exp(a(2) X(D))-Y(D; 
le 








置信 区 间 


nlparci: 非 线性 回归 参数 的 置信 区 间 
9] 革 人 |jel=iieilteisitc 卫 weekhzLES Te 
9 汪汪 1|erztielteisitc 表 Kijle 同 [=Uweleltc1ii 
ci=nlparci(.., alpha, alpha ) 

nlpredci: 非 线性 回归 预测 值 的 置信 区 间 
Wellsie 本 elsersieleitiiteleLsliigli 队 机 elsitc 古 [iie 隐 weieit2y 
Mete 本 el Eeeeldiiieieis ii 二 人 el 二 必 L=siie 风 EeelelU 
[.] = nlpredci(. "param1,val1, param2 ,val2,….) 

nlintool: 以 两 条 红线 绘制 预测 值 的 置信 区 间 
nlintool(x,y,fun,beta0) 





优化 函数 的 控制 





optimget: 得 到 最 优化 函数 选项 参数 值 
optimset: 创建 或 编辑 最 优化 函数 选项 参数 值 
调用 格 式 : 


opts=optimset(Cpara1 ,value1 ,para2 ,value2 ,….) 
将 指定 参数 para 赋 为 的 value 
jellliiSisl 
返回 所 有 可 以 供 optimset 设 置 的 参数 值 
el 二 eljelliiSa 
创建 一 个 结构 ， 所 有 的 值 都 被 赋予 [] 
jelKejeltiSiaitelieleljeik 尖 el LI 和 
从 oldopts 复 制 参数 值 ， 并 修改 相应 的 参数 值 
eljelk 二 eljelliSisit(ellelejelSaaIs ee 
oldqopts 和 newopts 组 合 ， 以 newopts 的 非 零 参 数值 履 盖 oldopts 中 的 对 应 值 





遗传 算法 与 直接 搜索 工具 箱 简介 


Matlab7.0 R14 以 后 的 版 本 中 包含 了 一 个 专门 设计 的 遗 
传 算 法 与 直接 搜索 工具 箱 (Genetic Algorithm and 
Direct Search Toolbox，GADST ) 


拓展 了 原 有 最 优化 工具 箱 的 功能 ， 可 以 求解 普通 方法 不 
钨 求解 的 问题 ， 如 查找 问题 


与 最 优化 工具 箱 紧密 结合 ， 可 以 充分 发 挥 工具 箱 功能 ， 
以 提高 求解 的 质量 。 对 于 某 些 问 题 可 以 获得 全 局 最 优 解 





遗传 算法 的 基本 概念 cp 


。 站 传 算法 是 一 种 模拟 自然 选择 、 生 物 进化 过 程 来 来 
解 问题 的 方法 
。 算法 要 点 
1 .创建 一 个 随机 种 群 (Populationm) 
2. 生 成 一 个 新 的 种 群 
al) 7 给 当前 种 群 的 每 一 个 个 体 
分 


b) 根据 适应 度 选择 父 厘 


C) 由 父 草 产生 子 昔 ， YL 交 
又 (Crossoven)、 变 异 (Mutationm) 和 迁移 (Migration) 进 行 


d) 用 子 草 蔡 换 当 前 种 群 ， 形 成 下 一 代 (Generation) 
3. 若 满足 停止 准则 ， 则 算法 停止 





。 在 命令 窗口 键入 gatool， 可 打开 


让 


了 ile Help 





IC 点 1LgDTItD 


i 卫 工 DD1 








了 itness function' 


Jamber of variables: 


Linear inequalities: 


Linear equalities: 


Bounds : 


Plot interval 





Best fitness 








_ | 了 xpectatiomn 





Seore diversity 








Stopping 











Custom function: 


| start 


Current gerneration: 


Status and results: 


世 


了 inal point: 





只 = b = 
上 Req 二 beq = 
Lower = Upper = 


onlinear constraint function' 






































LJjBest individual 了 ist 人 
_ | Genealogy 品 Range 
Seores 口 selec 

LNWax constraint 
Clear Stat 


Dptions:- 





| 


器 Population 





Fopulation type': 
Fopulation size: 


Creatiomn function: 


Initial population' 
Initial scores: 


Initial range: 


Double Vector 二 
20 

Uni form 二 
[] 

[] 

让 





由 了 itness scaling 








由 Selection 








由 Reproduction 








由 Iatation 








由 Crossoyer 








由 用 如 atiom 








由 委 gorithm settings 








国 ]ybrid funetion 








由 Stopping _ criteria 








由 Dutput _ function 








由 Display to _comman 


d_wimdow 














由 Yectorize 











Qui ck reference:- | < | 





Genetic Algorithm 


Solver 
Thistool corresponds to the ga 
function， 


Fitness function (required) is 
the objective function you wantto 
minimize.You can specify the 
function as afunction handle of 
the form @objfun, where 
objfun.misanv-flle that 
returns a SCalar. 


Number of variables (required) 
isthe number of independent 
variables forthe fitness function 


上 


上 


Options 
Specify options forthe Cenetic 
员 Igorithm solver. 


上 





RCI 





ee reE 本 eejliyuEeS 
f = exp(x(1))”(4*Xx(1)^A2+2*X(2)^A2+4*X(1)X(2)+2*X(2)+1); 


etelia [oa 区 e 放 切 丁 eeweiie4 
c1=[1.5+x(1)*Xx(2)-x(1)-Xx(2);-x(1T)*Xx(2)-10]; 
2 6 





函数 ga 的 用 法 (人 
调用 格式 : 


x=galfitnessfun, nvars, options) 











数学 模型 Xx=galfitnessfcn,nvars,A,b,Aed,beq,LB,UB,nonlcon;,options) 
[x, fval, reason, output, population, Scores] = gal(.…) 
输入 参数 : 
sssill 目标 函数 
nvars 解 向 量 的 长 度 
A,b 线性 不 等 式 约束 
ee 线性 等 式 约束 
LB，UB 自 变 量 的 下 、 上 限 
eliweli 非 线性 约束 
options gaoptimset 设 置 的 算法 参数 
输出 参数 : 
reason 表明 算法 终止 的 字符 
outpuit 表明 每 代 计 算 结果 和 算法 效果 的 结构 体 ， 


randstate: 随机 数 发 生 喜 的 状态 
randnstate: 均匀 分 别 随机 数 发 生 器 状态 
generations: 计算 代数 

funccount: 适应 度 函 数 的 计算 次 数 
message: 算法 的 终止 原因 

elejeil iiieli 每 列 是 最 终 的 种 群 
Scores 最 终 种 群 的 评价 





马 数 ga() 的 用 法 






TeilelE@IETXeisiiite72 


Xx0O=|[-1,1];a=[1,-1]j;b=1;a1=[1,1;b1=0; 

[x,fj=fmincon(C2objfun,x0,a,b,af1,b1, 赂 @nonlcom) 
eeEUeeiSistgaieiaei 肖 全 CieElelellelsSilCeEleleileissitilei 人 3 
CPULCSIeRelliel 人 是 eieje 下 :ieeigsj ecT(C2ele]jiigliE2 人 = 同 e 丙 0 罗 e 的 上 国 上 区 22201elal[ere1g 丰 oje1 
[x3,f3j=fmincon(Cobjfun,x2,a,b,a1,b1, 由 ,由 Ononlcon) 


Liel 站 elejiigliib 
f = exp(X(1))(4*X(1)A2+2*X(2)^A2+4”*X(1)X(2)+2*X(2)+1); 


function [c1,c2]=nonlcon(X) 
c1=[1.5+x(1)*x(2)-x(1)-Xx(2);-x(1)*Xx(2)-10]; 
C2=0; 


也 钨 、 1) 、 7 
误 传 算法 参数 设置 (全 
1. 图 形 参 数 

2. 

3. 适应 度 计 算 人 参数 -适应 度 测 量 值 变化 大 ， 则 高 测量 值 个 体 复制 快 ， 会 影响 


4. 
下 


6. 





种 群 会 数 -种 群 多 样 性 增加 ， 种 群 个 体 增加 可 能 提高 结果 质量 


搜索 空间 的 大 小 

选择 参数 

再 生 参 数 - 优 良 个体 多 则 可 以 使 最 适应 个 体 控制 种 群 ， 但 可 能 减少 搜索 有 
效 性 

变异 参数 -变异 和 交叉 是 遗传 算法 的 两 个 主要 过 程 ， 但 参数 的 设置 需 通过 
一 定 实 验 寻 找 最 佳 值 


7. 交叉 参数 
8. 
9. 


算法 设置 


10. 混 合 函 数 参 数 -~ 可 将 遗传 算法 与 其 它 算法 的 结合 ， 提 高 结果 
11. 停 止 条 件 参 数 

12. 输 出 函数 参数 

13. 显 示 到 命令 窗口 参数 

14. 向 量 参数 





直接 搜索 工具 箱 Ce 


。 与 着 传 算法 相同 ， 直 接 搜索 也 是 一 种 全 局 优化 算法 


”在 命令 窗口 键入 psearchtool， 可 打开 模式 搜索 GUI 
求解 问题 ， 其 使 用 方法 与 遗传 算法 工具 箱 类 似 


。 直接 搜索 算法 的 核心 函数 是 patternsearch 





。 采 用 图 表 直 观 观察 数据 是 否 有 异常 ， 对 于 异常 数 
据点 应 舍弃 ， 人 必要 时 应 重复 实验 

。 数 据 是 否 一 臻 

。 数 据 是 否 不 足 

"应 重视 实验 设计 





二 2 \ 多 人 一 三 、 人 入 让 
间 区 反应 器 中 串联 -平行 复杂 反应 系统 


在 间歇 反应 器 中 进行 液 相 反应 制备 产物 B， 反 应 网 络 如 图 5-1 所 示 。 反 
应 可 在 180~260C 的 温度 范围 内 进行 ， 反 应 物 X 大 量 过 剩 ， 而 C, D 和 
E 为 副 产 物 。 各 反应 均 为 一 级 动力 学 关系 : Tf= - kC， 式 中 辐 喝 


已 知 参 数 : k01 = 5.78052 x 1010，k02 = 3.92317 x 1012，k03 = 

1.64254 x 104，k04 = 6.264 x 108，Ea1=124670，Ea2=150386， 
Ea3=77954，Ea4=111528。 初 始 浓度 : CA=1kmoym3， 其 余 物 质 浓度 为 0。 
用 最 优化 方法 计算 使 产物 B 收 率 最 大 的 最 优 反 应 温度 。 











间 歌 反应 器 中 串联 -平行 复杂 反应 系统 








weliEeit2VYXeisliieis 
R = 8.31434; ELSEweliSitcIl 人 Te 
= [5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8]; 
= [124670 150386 77954 111528]; % Activation energy, kJ/kmol 
=[10000; os%lnitialConcentration,C0(), kmoMmA^3 
tspan = [0 1e4]; 
T0O = 200+2793; % Reactor temperature, Kelvin 
T=fminsearch(@OOobjFunc,T0,[,tspan,C0O,k0,Ea,R); 
fprintf(\t 使 产物 B 收 率 最 大 的 最 优 反应 温度 为 %3.1fCNn',T - 273) 
26 -一 
function f = ObjFunc(T,tspan,C0,k0,Ea,R) 
[t,C] = ode45(@MassEquations, tspan, CO0,[,k0,Ea,R,T); 
CBmax = max(C(:,2)); 
f= -CBmax: 


function dCdt = MassEquations(t,C,k0,Ea,R,T) 
k = k0.*exp(-Ea/(R*T)); 

k(5) = 2.16667E-04; 

rA = -(k(1)+K(2))C(T); 

rB = k(1T) 六 CD)-k(3 六 C 
rC = k(2) 六 CD)-K(4) 
rD = K(3) CC(2)-K(5) 
rE = K(4) 六 CC(3)+K(5) 六 
dcCdt = [TA; rB; rC;r 


疝 寺 


WE 





乙烯 加 气 反 应 动力 学 参数 估计 
乙烯 (E) 加 气 (H) 制 乙 烷 (EA) 反应 如 下 上 事 了 2 玉 风 2 月 


微分 反应 器 中 测 得 的 实验 数据 如 下 : 反应 速率 方程 为 ; 


palatm 








om | 根据 吉 中 效 据 痛 定 


动力 学 参数 Kk 和 KF 





EECZNRETNETT 有 IETT 
IIETNETTNETT 有 


一 
一 
7 





aie eeLsiiie7 

rA = [1.04 3.13 5.21 3.82 4.19 2.391 3.867 2.199 0.75]; 
PRe 二 | 中 省 8505050505|: 

民 几 宇 几 Eee 

9%6 线性 拟 合 

RA = Pe.*Ph.rA; y = RA':;X = [ones(size(y)) Pe ]; 

b = X\VY;k= 1/b(1); 

KE = b(2)kK; 

% 非 线性 拟 合 

1 下 


eaERL=SiiLeldesilelliIRs>ditiEUeEelilielila[Uleie[ 百 =Uwelej[1]E 二 汪汪 


Isqnonlin(@ObjFunc,beta0, 几 名 , 册 ,rA,Pe,Ph) 
有 er=Tielileisit 机 [sjlelil EUweleiti 
9% 残 差 关于 拟 合 值 的 残 差 图 
rAc = RateA(beta,Pe,Ph) 
el eltgraealsrsile[ 引 -| 
xlabel( 反 应 速率 ( -frA) 拟 合 值 , molW(kg cat. S)) 
ylabel(' 残 差 R, mol(kg cat. s)) 
refline(0,0) 
9% 参数 辨识 结果 
fprintf(CEstimated Parameters An') 
fprintf(Ntk = %.3f iA %.3fwm',beta(1),ci(1,2)-beta(1)) 
fprintf(NKE = %.3f iA %.3fwm',beta(2),ci(2,2)-beta(2)) 
fprintf(AtThe sum of the squares is: %.3f ,resnorm) 
站 
function f= ObjFunc(beta,rA,Pe,Ph) 
f=rA-RateAlbeta,Pe,Ph); 
本 
function rA = RateA(beta,Pe,Ph) 
rA = betal(1)*Pe.Ph./(1+beta(2)Pe); 


祝 大 家 在 今后 的 学 习 生 活 中 一 切 顺 利 ! 


